1 |
chrisfen |
1287 |
/*************************************************************************** |
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 |
|
|
*/ |