1 |
/*************************************************************************** |
2 |
************************************************************************** |
3 |
|
4 |
S2kit 1.0 |
5 |
|
6 |
A lite version of Spherical Harmonic Transform Kit |
7 |
|
8 |
Peter Kostelec, Dan Rockmore |
9 |
{geelong,rockmore}@cs.dartmouth.edu |
10 |
|
11 |
Contact: Peter Kostelec |
12 |
geelong@cs.dartmouth.edu |
13 |
|
14 |
Copyright 2004 Peter Kostelec, Dan Rockmore |
15 |
|
16 |
This file is part of S2kit. |
17 |
|
18 |
S2kit is free software; you can redistribute it and/or modify |
19 |
it under the terms of the GNU General Public License as published by |
20 |
the Free Software Foundation; either version 2 of the License, or |
21 |
(at your option) any later version. |
22 |
|
23 |
S2kit is distributed in the hope that it will be useful, |
24 |
but WITHOUT ANY WARRANTY; without even the implied warranty of |
25 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
26 |
GNU General Public License for more details. |
27 |
|
28 |
You should have received a copy of the GNU General Public License |
29 |
along with S2kit; if not, write to the Free Software |
30 |
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
31 |
|
32 |
See the accompanying LICENSE file for details. |
33 |
|
34 |
************************************************************************ |
35 |
************************************************************************/ |
36 |
|
37 |
|
38 |
/* |
39 |
source file which contains the function that generates the |
40 |
weights for a bandwidth bw legendre transform. Basically, |
41 |
it contains the implementation of the formula as defined in |
42 |
the tensor paper, and also given in the so(3) paper. It's |
43 |
just mentioned in the s^2 paper! |
44 |
|
45 |
This formula is slightly different from the one given in the |
46 |
original DH paper because they were sampling at the poles in |
47 |
that paper, and now we're not. |
48 |
|
49 |
In pseudo-TeX, the formula for the bandwidth B weights is |
50 |
|
51 |
w_B(j) = 2/B sin((pi*(2j+1))/(4B)) * |
52 |
sum_{k=0}^{B-1} 1/(2k+1)*sin((2j+1)(2k+1)pi/(4B)) |
53 |
|
54 |
where j = 0, 1, ..., 2 B - 1 |
55 |
|
56 |
|
57 |
Note that if you want to use these weights for an *odd* |
58 |
order transform, given the way the code is set up, you |
59 |
have to MULTIPLY the j-th weight by sin(pi*(2j+1)/(4B)) |
60 |
|
61 |
*/ |
62 |
|
63 |
#include <math.h> |
64 |
|
65 |
/* |
66 |
makeweights: given a bandwidth bw, make weights for |
67 |
both even *and* odd order Legendre transforms. |
68 |
|
69 |
bw -> bandwidth of transform |
70 |
weights -> pointer to double array of length 4*bw; this |
71 |
array will contain the even and odd weights; |
72 |
even weights start at weight[0], and odd weights |
73 |
start at weights[2*bw] |
74 |
|
75 |
*/ |
76 |
|
77 |
void makeweights( int bw, |
78 |
double *weights ) |
79 |
{ |
80 |
int j, k ; |
81 |
double fudge ; |
82 |
double tmpsum ; |
83 |
|
84 |
fudge = M_PI/((double)(4*bw)) ; |
85 |
|
86 |
|
87 |
for ( j = 0 ; j < 2*bw ; j ++ ) |
88 |
{ |
89 |
tmpsum = 0.0 ; |
90 |
for ( k = 0 ; k < bw ; k ++ ) |
91 |
tmpsum += 1./((double)(2*k+1)) * |
92 |
sin((double)((2*j+1)*(2*k+1))*fudge); |
93 |
tmpsum *= sin((double)(2*j+1)*fudge); |
94 |
tmpsum *= 2./((double) bw) ; |
95 |
|
96 |
weights[j] = tmpsum ; |
97 |
weights[j + 2*bw] = tmpsum * sin((double)(2*j+1)*fudge); |
98 |
} |
99 |
|
100 |
} |
101 |
|
102 |
/********************************************************/ |
103 |
/********************************************************/ |
104 |
/********************************************************/ |
105 |
/********************************************************/ |
106 |
|
107 |
/* just a hack to test the above function */ |
108 |
/* |
109 |
|
110 |
#include <stdio.h> |
111 |
#include <stdlib.h> |
112 |
|
113 |
int main( int argc, |
114 |
char **argv ) |
115 |
{ |
116 |
int i, bw ; |
117 |
double *weights ; |
118 |
|
119 |
|
120 |
bw = atoi(argv[1]); |
121 |
|
122 |
weights = (double *)malloc(sizeof(double)*4*bw); |
123 |
|
124 |
makeweights(bw, weights); |
125 |
|
126 |
for ( i = 0 ; i < 2*bw ; i ++ ) |
127 |
printf("%d\t%f\t%f\n", |
128 |
i, weights[i], weights[i+2*bw]); |
129 |
|
130 |
free( weights ); |
131 |
|
132 |
return 0 ; |
133 |
|
134 |
} |
135 |
|
136 |
*/ |