codebase / stimlib

Commit eaf2cb076d7dfca096a820c1530dbad59872366b

Authored by David Mayerich
1 parent ba72c5b0

renamed spherical harmonics header

Showing 1 changed file with 236 additions and 236 deletions   Show diff stats
Inline Side-by-side
stim/visualization/sph_harmonics.h renamed to stim/visualization/spharmonics.h
  Show/Hide comments   View file @eaf2cb0
1   -#ifndef STIM_SPH_HARMONICS
2   -#define STIM_SPH_HARMONICS
3   -
4   -#include <GL/gl.h>
5   -
6   -#include <stim/gl/error.h>
7   -#include <stim/visualization/colormap.h>
8   -#include <vector>
9   -
10   -#define PI 3.14159
11   -#define WIRE_SCALE 1.001
12   -namespace stim{
13   -
14   - class sph_harmonics{
15   -
16   - private:
17   -
18   - double* func; //stores the raw function data (samples at each point)
19   -
  1 +#ifndef STIM_SPH_HARMONICS
  2 +#define STIM_SPH_HARMONICS
  3 +
  4 +#include <GL/gl.h>
  5 +
  6 +#include <stim/gl/error.h>
  7 +#include <stim/visualization/colormap.h>
  8 +#include <vector>
  9 +
  10 +#define PI 3.14159
  11 +#define WIRE_SCALE 1.001
  12 +namespace stim{
  13 +
  14 + class sph_harmonics{
  15 +
  16 + private:
  17 +
  18 + double* func; //stores the raw function data (samples at each point)
  19 +
20 20 GLuint color_tex; //texture map that acts as a colormap for the spherical function
21 21  
22   - unsigned int N; //resolution of the spherical grid
23   -
24   - std::vector<double> C; //list of SH coefficients
25   -
26   -
27   - //evaluates an associated Legendre polynomial (-l <= m <= l)
28   - double P(int l,int m,double x)
29   - {
30   - // evaluate an Associated Legendre Polynomial P(l,m,x) at x
31   - double pmm = 1.0;
32   - if(m>0) {
33   - double somx2 = sqrt((1.0-x)*(1.0+x));
34   - double fact = 1.0;
35   - for(int i=1; i<=m; i++) {
36   - pmm *= (-fact) * somx2;
37   - fact += 2.0;
38   - }
39   - }
40   - if(l==m) return pmm;
41   - double pmmp1 = x * (2.0*m+1.0) * pmm;
42   - if(l==m+1) return pmmp1;
43   - double pll = 0.0;
44   - for(int ll=m+2; ll<=l; ++ll) {
45   - pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
46   - pmm = pmmp1;
47   - pmmp1 = pll;
48   - }
49   - return pll;
50   - }
51   -
52   - //recursively calculate a factorial given a positive integer n
53   - unsigned int factorial(unsigned int n) {
54   - if (n == 0)
55   - return 1;
56   - return n * factorial(n - 1);
57   - }
58   -
59   - //calculate the SH scaling constant
60   - double K(int l, int m){
61   -
62   - // renormalisation constant for SH function
63   - double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
64   - return sqrt(temp);
65   - }
66   -
67   - //calculate the value of the SH basis function (l, m) at (theta, phi)
68   - //here, theta = [0, PI], phi = [0, 2*PI]
69   - double SH(int l, int m, double theta, double phi){
70   - // return a point sample of a Spherical Harmonic basis function
71   - // l is the band, range [0..N]
72   - // m in the range [-l..l]
73   - // theta in the range [0..Pi]
74   - // phi in the range [0..2*Pi]
75   - const double sqrt2 = sqrt(2.0);
76   - if(m==0) return K(l,0)*P(l,m,cos(theta));
77   - else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
78   - else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
79   - }
80   -
81   - void gen_function(){
82   -
83   - //initialize the function to zero
84   - memset(func, 0, sizeof(double) * N * N);
85   -
86   - double theta, phi;
87   - double result;
88   - int l, m;
89   -
90   - l = m = 0;
91   - for(unsigned int c = 0; c < C.size(); c++){
92   -
93   -
  22 + unsigned int N; //resolution of the spherical grid
  23 +
  24 + std::vector<double> C; //list of SH coefficients
  25 +
  26 +
  27 + //evaluates an associated Legendre polynomial (-l <= m <= l)
  28 + double P(int l,int m,double x)
  29 + {
  30 + // evaluate an Associated Legendre Polynomial P(l,m,x) at x
  31 + double pmm = 1.0;
  32 + if(m>0) {
  33 + double somx2 = sqrt((1.0-x)*(1.0+x));
  34 + double fact = 1.0;
  35 + for(int i=1; i<=m; i++) {
  36 + pmm *= (-fact) * somx2;
  37 + fact += 2.0;
  38 + }
  39 + }
  40 + if(l==m) return pmm;
  41 + double pmmp1 = x * (2.0*m+1.0) * pmm;
  42 + if(l==m+1) return pmmp1;
  43 + double pll = 0.0;
  44 + for(int ll=m+2; ll<=l; ++ll) {
  45 + pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
  46 + pmm = pmmp1;
  47 + pmmp1 = pll;
  48 + }
  49 + return pll;
  50 + }
  51 +
  52 + //recursively calculate a factorial given a positive integer n
  53 + unsigned int factorial(unsigned int n) {
  54 + if (n == 0)
  55 + return 1;
  56 + return n * factorial(n - 1);
  57 + }
  58 +
  59 + //calculate the SH scaling constant
  60 + double K(int l, int m){
  61 +
  62 + // renormalisation constant for SH function
  63 + double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
  64 + return sqrt(temp);
  65 + }
  66 +
  67 + //calculate the value of the SH basis function (l, m) at (theta, phi)
  68 + //here, theta = [0, PI], phi = [0, 2*PI]
  69 + double SH(int l, int m, double theta, double phi){
  70 + // return a point sample of a Spherical Harmonic basis function
  71 + // l is the band, range [0..N]
  72 + // m in the range [-l..l]
  73 + // theta in the range [0..Pi]
  74 + // phi in the range [0..2*Pi]
  75 + const double sqrt2 = sqrt(2.0);
  76 + if(m==0) return K(l,0)*P(l,m,cos(theta));
  77 + else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
  78 + else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
  79 + }
  80 +
  81 + void gen_function(){
  82 +
  83 + //initialize the function to zero
  84 + memset(func, 0, sizeof(double) * N * N);
  85 +
  86 + double theta, phi;
  87 + double result;
  88 + int l, m;
  89 +
  90 + l = m = 0;
  91 + for(unsigned int c = 0; c < C.size(); c++){
  92 +
  93 +
94 94 for(unsigned int xi = 0; xi < N; xi++)
95 95 for(unsigned int yi = 0; yi < N; yi++){
96 96  
... ... @@ -98,28 +98,28 @@ namespace stim{
98 98 phi = PI * ((double)yi / (N-1));
99 99 result = C[c] * SH(l, m, phi, theta); //phi and theta are reversed here (damn physicists)
100 100 func[yi * N + xi] += result;
101   - }
102   -
103   - m++; //increment m
104   -
105   - //if we're in a new tier, increment l and set m = -l
106   - if(m > l){
107   - l++;
108   - m = -l;
109   - }
110   - }
111   - }
112   -
113   - void gl_prep_draw(){
114   -
115   - //enable depth testing
116   - //this has to be used instead of culling because the sphere can have negative values
117   - glEnable(GL_DEPTH_TEST);
118   - glDepthMask(GL_TRUE);
119   - glEnable(GL_TEXTURE_2D); //enable 2D texture mapping
120   - }
121   -
122   - //draw a texture mapped sphere representing the function surface
  101 + }
  102 +
  103 + m++; //increment m
  104 +
  105 + //if we're in a new tier, increment l and set m = -l
  106 + if(m > l){
  107 + l++;
  108 + m = -l;
  109 + }
  110 + }
  111 + }
  112 +
  113 + void gl_prep_draw(){
  114 +
  115 + //enable depth testing
  116 + //this has to be used instead of culling because the sphere can have negative values
  117 + glEnable(GL_DEPTH_TEST);
  118 + glDepthMask(GL_TRUE);
  119 + glEnable(GL_TEXTURE_2D); //enable 2D texture mapping
  120 + }
  121 +
  122 + //draw a texture mapped sphere representing the function surface
123 123 void gl_draw_sphere() {
124 124  
125 125 //PI is used to convert from spherical to cartesian coordinates
... ... @@ -169,133 +169,133 @@ namespace stim{
169 169 }
170 170 glEnd();
171 171 }
172   - }
173   -
174   - //draw a wire frame sphere representing the function surface
175   - void gl_draw_wireframe() {
176   -
177   - //PI is used to convert from spherical to cartesian coordinates
178   - //const double PI = 3.14159;
179   -
180   - //bind the 2D texture representing the color map
181   - glDisable(GL_TEXTURE_2D);
182   - glColor3f(0.0f, 0.0f, 0.0f);
183   -
184   - //Draw the Sphere
185   - int i, j;
186   -
187   - for(i = 1; i <= N-1; i++) {
188   - double phi0 = PI * ((double) (i - 1) / (N-1));
189   - double phi1 = PI * ((double) i / (N-1));
190   -
191   - glBegin(GL_LINE_STRIP);
192   - for(j = 0; j <= N; j++) {
193   -
194   - //calculate the indices into the function array
195   - int phi0_i = i-1;
196   - int phi1_i = i;
197   - int theta_i = j;
198   - if(theta_i == N)
199   - theta_i = 0;
200   -
201   - double v0 = func[phi0_i * N + theta_i];
202   - double v1 = func[phi1_i * N + theta_i];
203   -
204   - v0 = fabs(v0);
205   - v1 = fabs(v1);
206   -
207   -
208   - double theta = 2 * PI * (double) (j - 1) / N;
209   - double x0 = WIRE_SCALE * v0 * cos(theta) * sin(phi0);
210   - double y0 = WIRE_SCALE * v0 * sin(theta) * sin(phi0);
211   - double z0 = WIRE_SCALE * v0 * cos(phi0);
212   -
213   - double x1 = WIRE_SCALE * v1 * cos(theta) * sin(phi1);
214   - double y1 = WIRE_SCALE * v1 * sin(theta) * sin(phi1);
215   - double z1 = WIRE_SCALE * v1 * cos(phi1);
216   -
217   - glTexCoord2f(theta / (2 * PI), phi0 / PI);
218   - glVertex3f(x0, y0, z0);
219   -
220   - glTexCoord2f(theta / (2 * PI), phi1 / PI);
221   - glVertex3f(x1, y1, z1);
222   - }
223   - glEnd();
224   - }
225   - }
226   -
227   - void init(unsigned int n){
228   -
229   - //set the sphere resolution
230   - N = n;
231   -
232   - //allocate space for the color map
233   - unsigned int bytes = N * N * sizeof(unsigned char) * 3;
  172 + }
  173 +
  174 + //draw a wire frame sphere representing the function surface
  175 + void gl_draw_wireframe() {
  176 +
  177 + //PI is used to convert from spherical to cartesian coordinates
  178 + //const double PI = 3.14159;
  179 +
  180 + //bind the 2D texture representing the color map
  181 + glDisable(GL_TEXTURE_2D);
  182 + glColor3f(0.0f, 0.0f, 0.0f);
  183 +
  184 + //Draw the Sphere
  185 + int i, j;
  186 +
  187 + for(i = 1; i <= N-1; i++) {
  188 + double phi0 = PI * ((double) (i - 1) / (N-1));
  189 + double phi1 = PI * ((double) i / (N-1));
  190 +
  191 + glBegin(GL_LINE_STRIP);
  192 + for(j = 0; j <= N; j++) {
  193 +
  194 + //calculate the indices into the function array
  195 + int phi0_i = i-1;
  196 + int phi1_i = i;
  197 + int theta_i = j;
  198 + if(theta_i == N)
  199 + theta_i = 0;
  200 +
  201 + double v0 = func[phi0_i * N + theta_i];
  202 + double v1 = func[phi1_i * N + theta_i];
  203 +
  204 + v0 = fabs(v0);
  205 + v1 = fabs(v1);
  206 +
  207 +
  208 + double theta = 2 * PI * (double) (j - 1) / N;
  209 + double x0 = WIRE_SCALE * v0 * cos(theta) * sin(phi0);
  210 + double y0 = WIRE_SCALE * v0 * sin(theta) * sin(phi0);
  211 + double z0 = WIRE_SCALE * v0 * cos(phi0);
  212 +
  213 + double x1 = WIRE_SCALE * v1 * cos(theta) * sin(phi1);
  214 + double y1 = WIRE_SCALE * v1 * sin(theta) * sin(phi1);
  215 + double z1 = WIRE_SCALE * v1 * cos(phi1);
  216 +
  217 + glTexCoord2f(theta / (2 * PI), phi0 / PI);
  218 + glVertex3f(x0, y0, z0);
  219 +
  220 + glTexCoord2f(theta / (2 * PI), phi1 / PI);
  221 + glVertex3f(x1, y1, z1);
  222 + }
  223 + glEnd();
  224 + }
  225 + }
  226 +
  227 + void init(unsigned int n){
  228 +
  229 + //set the sphere resolution
  230 + N = n;
  231 +
  232 + //allocate space for the color map
  233 + unsigned int bytes = N * N * sizeof(unsigned char) * 3;
234 234 unsigned char* color_image;
235   - color_image = (unsigned char*) malloc(bytes);
236   -
237   - //allocate space for the function
238   - func = (double*) malloc(N * N * sizeof(double));
239   -
240   - //generate a function (temporary)
241   - gen_function();
242   -
243   - //generate a colormap from the function
244   - stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);
245   -
246   - //prep everything for drawing
247   - gl_prep_draw();
248   -
249   - //generate an OpenGL texture map in the current context
  235 + color_image = (unsigned char*) malloc(bytes);
  236 +
  237 + //allocate space for the function
  238 + func = (double*) malloc(N * N * sizeof(double));
  239 +
  240 + //generate a function (temporary)
  241 + gen_function();
  242 +
  243 + //generate a colormap from the function
  244 + stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);
  245 +
  246 + //prep everything for drawing
  247 + gl_prep_draw();
  248 +
  249 + //generate an OpenGL texture map in the current context
250 250 glGenTextures(1, &color_tex);
251 251 //bind the texture
252   - glBindTexture(GL_TEXTURE_2D, color_tex);
253   -
254   - //copy the color data from the buffer to the GPU
255   - glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
256   -
257   - //initialize all of the texture parameters
  252 + glBindTexture(GL_TEXTURE_2D, color_tex);
  253 +
  254 + //copy the color data from the buffer to the GPU
  255 + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
  256 +
  257 + //initialize all of the texture parameters
258 258 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
259 259 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
260 260 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
261 261 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
262   - glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
263   -
264   - //free the buffer
265   - free(color_image);
266   - }
267   -
268   -
269   - public:
270   -
271   - void glRender(){
272   - //set all OpenGL parameters required for drawing
273   - gl_prep_draw();
274   -
275   - //draw the sphere
276   - gl_draw_sphere();
277   - //gl_draw_wireframe();
278   -
279   - }
280   -
281   - void glInit(unsigned int n){
282   - init(n);
283   - }
284   -
285   - void push(double c){
286   - C.push_back(c);
287   - }
288   -
289   -
290   -
291   -
292   -
293   - }; //end class sph_harmonics
294   -
295   -
296   -
297   -
298   -}
299   -
300   -
301   -#endif
  262 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
  263 +
  264 + //free the buffer
  265 + free(color_image);
  266 + }
  267 +
  268 +
  269 + public:
  270 +
  271 + void glRender(){
  272 + //set all OpenGL parameters required for drawing
  273 + gl_prep_draw();
  274 +
  275 + //draw the sphere
  276 + gl_draw_sphere();
  277 + //gl_draw_wireframe();
  278 +
  279 + }
  280 +
  281 + void glInit(unsigned int n){
  282 + init(n);
  283 + }
  284 +
  285 + void push(double c){
  286 + C.push_back(c);
  287 + }
  288 +
  289 +
  290 +
  291 +
  292 +
  293 + }; //end class sph_harmonics
  294 +
  295 +
  296 +
  297 +
  298 +}
  299 +
  300 +
  301 +#endif
... ...