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 |
|
|
/******************************************************************** |
40 |
|
|
|
41 |
|
|
FST_semi_memo.c - routines to perform convolutions on the |
42 |
|
|
2-sphere using a combination of semi-naive and naive algorithms. |
43 |
|
|
|
44 |
|
|
ASSUMES THAT ALL PRECOMPUTED-DATA IS IN MEMORY, AND NOT TO BE |
45 |
|
|
READ FROM THE DISK. |
46 |
|
|
|
47 |
|
|
The primary functions in this package are |
48 |
|
|
|
49 |
|
|
1) FST_semi_memo() - computes the spherical harmonic expansion. |
50 |
|
|
2) InvFST_semi_memo() - computes the inverse spherical harmonic transform. |
51 |
|
|
3) FZT_semi_memo() - computes the zonal harmonic transform. |
52 |
|
|
4) TransMult() - Multiplies harmonic coefficients using Driscoll-Healy |
53 |
|
|
result. Dual of convolution in "time" domain. |
54 |
|
|
5) Conv2Sphere_semi_memo() - Convolves two functins defined on the 2-sphere, |
55 |
|
|
using seminaive transform |
56 |
|
|
|
57 |
|
|
and one utility function: |
58 |
|
|
|
59 |
|
|
1) seanindex(): Given bandwidth bw, seanindex(m,l,bw) will give the |
60 |
|
|
position of the coefficient f-hat(m,l) in the one-row array |
61 |
|
|
|
62 |
|
|
For descriptions on calling these functions, see the documentation |
63 |
|
|
preceding each function. |
64 |
|
|
|
65 |
|
|
*/ |
66 |
|
|
|
67 |
|
|
#include <math.h> |
68 |
|
|
#include <stdio.h> |
69 |
|
|
#include <stdlib.h> |
70 |
|
|
#include <string.h> |
71 |
|
|
|
72 |
|
|
#include "fftw3.h" |
73 |
|
|
|
74 |
|
|
#include "makeweights.h" |
75 |
|
|
#include "cospmls.h" |
76 |
|
|
#include "primitive.h" |
77 |
|
|
#include "naive_synthesis.h" |
78 |
|
|
#include "seminaive.h" |
79 |
|
|
#include "FST_semi_memo.h" |
80 |
|
|
|
81 |
|
|
|
82 |
|
|
|
83 |
|
|
/************************************************************************/ |
84 |
|
|
|
85 |
|
|
|
86 |
|
|
/***************************************************************** |
87 |
|
|
|
88 |
|
|
Given bandwidth bw, seanindex(m,l,bw) will give the position of the |
89 |
|
|
coefficient f-hat(m,l) in the one-row array that Sean stores the spherical |
90 |
|
|
coefficients. This is needed to help preserve the symmetry that the |
91 |
|
|
coefficients have: (l = degree, m = order, and abs(m) <= l) |
92 |
|
|
|
93 |
|
|
f-hat(l,-m) = (-1)^m * conjugate( f-hat(l,m) ) |
94 |
|
|
|
95 |
|
|
Thanks for your help Mark! |
96 |
|
|
|
97 |
|
|
******************************************************************/ |
98 |
|
|
|
99 |
|
|
|
100 |
|
|
int seanindex(int m, |
101 |
|
|
int l, |
102 |
|
|
int bw) |
103 |
|
|
{ |
104 |
|
|
int bigL; |
105 |
|
|
|
106 |
|
|
bigL = bw - 1; |
107 |
|
|
|
108 |
|
|
if( m >= 0 ) |
109 |
|
|
return( m * ( bigL + 1 ) - ( ( m * (m - 1) ) /2 ) + ( l - m ) ); |
110 |
|
|
else |
111 |
|
|
return( ( ( bigL * ( bigL + 3 ) ) /2 ) + 1 + |
112 |
|
|
( ( bigL + m ) * ( bigL + m + 1 ) / 2 ) + ( l - abs( m ) ) ); |
113 |
|
|
} |
114 |
|
|
|
115 |
|
|
|
116 |
|
|
/************************************************************************/ |
117 |
|
|
/* performs a spherical harmonic transform using the semi-naive |
118 |
|
|
and naive algorithms |
119 |
|
|
|
120 |
|
|
bw -> bandwidth of problem |
121 |
|
|
size -> size = 2*bw -> dimension of input array (recall that |
122 |
|
|
sampling is done at twice the bandwidth) |
123 |
|
|
|
124 |
|
|
The inputs rdata and idata are expected to be pointers to |
125 |
|
|
size x size arrays. The array rdata contains the real parts |
126 |
|
|
of the function samples, and idata contains the imaginary |
127 |
|
|
parts. |
128 |
|
|
|
129 |
|
|
rcoeffs and icoeffs are expected to be pointers to bw x bw arrays, |
130 |
|
|
and will contain the harmonic coefficients in a "linearized" form. |
131 |
|
|
The array rcoeffs contains the real parts of the coefficients, |
132 |
|
|
and icoeffs contains the imaginary parts. |
133 |
|
|
|
134 |
|
|
spharmonic_pml_table should be a (double **) pointer to |
135 |
|
|
the result of a call to Spharmonic_Pml_Table. Because this |
136 |
|
|
table is re-used in the inverse transform, and because for |
137 |
|
|
timing purposes the computation of the table is not included, |
138 |
|
|
it is passed in as an argument. Also, at some point this |
139 |
|
|
code may be used as par of a series of convolutions, so |
140 |
|
|
reducing repetitive computation is prioritized. |
141 |
|
|
|
142 |
|
|
spharmonic_pml_table will be an array of (double *) pointers |
143 |
|
|
the array being of length TableSize(m,bw) |
144 |
|
|
|
145 |
|
|
workspace needs to be a double pointer to an array of size |
146 |
|
|
(8 * bw^2) + (7 * bw). |
147 |
|
|
|
148 |
|
|
cutoff -> what order to switch from semi-naive to naive |
149 |
|
|
algorithm. |
150 |
|
|
|
151 |
|
|
dataformat =0 -> samples are complex, =1 -> samples real |
152 |
|
|
|
153 |
|
|
|
154 |
|
|
Output Ordering of coeffs f(m,l) is |
155 |
|
|
f(0,0) f(0,1) f(0,2) ... f(0,bw-1) |
156 |
|
|
f(1,1) f(1,2) ... f(1,bw-1) |
157 |
|
|
etc. |
158 |
|
|
f(bw-2,bw-2), f(bw-2,bw-1) |
159 |
|
|
f(bw-1,bw-1) |
160 |
|
|
f(-(bw-1),bw-1) |
161 |
|
|
f(-(bw-2),bw-2) f(-(bw-2),bw-1) |
162 |
|
|
etc. |
163 |
|
|
f(-2,2) ... f(-2,bw-1) |
164 |
|
|
f(-1,1) f(-1,2) ... f(-1,bw-1) |
165 |
|
|
|
166 |
|
|
|
167 |
|
|
This only requires an array of size (bw*bw). If zero-padding |
168 |
|
|
is used to make the indexing nice, then you need a an |
169 |
|
|
(2bw-1) * bw array - but that is not done here. |
170 |
|
|
Because of the amount of space necessary for doing |
171 |
|
|
large transforms, it is important not to use any |
172 |
|
|
more than necessary. |
173 |
|
|
|
174 |
|
|
*/ |
175 |
|
|
|
176 |
|
|
void FST_semi_memo(double *rdata, double *idata, |
177 |
|
|
double *rcoeffs, double *icoeffs, |
178 |
|
|
int bw, |
179 |
|
|
double **seminaive_naive_table, |
180 |
|
|
double *workspace, |
181 |
|
|
int dataformat, |
182 |
|
|
int cutoff, |
183 |
|
|
fftw_plan *dctPlan, |
184 |
|
|
fftw_plan *fftPlan, |
185 |
|
|
double *weights ) |
186 |
|
|
{ |
187 |
|
|
int size, m, i, j; |
188 |
|
|
int dummy0, dummy1 ; |
189 |
|
|
double *rres, *ires; |
190 |
|
|
double *rdataptr, *idataptr; |
191 |
|
|
double *fltres, *scratchpad; |
192 |
|
|
double *eval_pts; |
193 |
|
|
double pow_one; |
194 |
|
|
double tmpA, tmpSize ; |
195 |
|
|
|
196 |
|
|
size = 2*bw ; |
197 |
|
|
tmpSize = 1./ ((double ) size); |
198 |
|
|
tmpA = sqrt( 2. * M_PI ) ; |
199 |
|
|
|
200 |
|
|
rres = workspace; /* needs (size * size) = (4 * bw^2) */ |
201 |
|
|
ires = rres + (size * size); /* needs (size * size) = (4 * bw^2) */ |
202 |
|
|
fltres = ires + (size * size) ; /* needs bw */ |
203 |
|
|
eval_pts = fltres + bw; /* needs (2*bw) */ |
204 |
|
|
scratchpad = eval_pts + (2*bw); /* needs (4 * bw) */ |
205 |
|
|
|
206 |
|
|
/* total workspace is (8 * bw^2) + (7 * bw) */ |
207 |
|
|
|
208 |
|
|
/* do the FFTs along phi */ |
209 |
|
|
fftw_execute_split_dft( *fftPlan, |
210 |
|
|
rdata, idata, |
211 |
|
|
rres, ires ) ; |
212 |
|
|
/* |
213 |
|
|
normalize |
214 |
|
|
|
215 |
|
|
note that I'm getting the sqrt(2pi) in there at |
216 |
|
|
this point ... to account for the fact that the spherical |
217 |
|
|
harmonics are of norm 1: I need to account for |
218 |
|
|
the fact that the associated Legendres are |
219 |
|
|
of norm 1 |
220 |
|
|
*/ |
221 |
|
|
tmpSize *= tmpA ; |
222 |
|
|
for( j = 0 ; j < size*size ; j ++ ) |
223 |
|
|
{ |
224 |
|
|
rres[j] *= tmpSize ; |
225 |
|
|
ires[j] *= tmpSize ; |
226 |
|
|
} |
227 |
|
|
|
228 |
|
|
/* point to start of output data buffers */ |
229 |
|
|
rdataptr = rcoeffs; |
230 |
|
|
idataptr = icoeffs; |
231 |
|
|
|
232 |
|
|
for (m=0; m<bw; m++) { |
233 |
|
|
/* |
234 |
|
|
fprintf(stderr,"m = %d\n",m); |
235 |
|
|
*/ |
236 |
|
|
|
237 |
|
|
/*** test to see if before cutoff or after ***/ |
238 |
|
|
if (m < cutoff){ |
239 |
|
|
|
240 |
|
|
/* do the real part */ |
241 |
|
|
SemiNaiveReduced(rres+(m*size), |
242 |
|
|
bw, |
243 |
|
|
m, |
244 |
|
|
fltres, |
245 |
|
|
scratchpad, |
246 |
|
|
seminaive_naive_table[m], |
247 |
|
|
weights, |
248 |
|
|
dctPlan); |
249 |
|
|
|
250 |
|
|
/* now load real part of coefficients into output space */ |
251 |
|
|
memcpy(rdataptr, fltres, sizeof(double) * (bw - m)); |
252 |
|
|
|
253 |
|
|
rdataptr += bw-m; |
254 |
|
|
|
255 |
|
|
/* do imaginary part */ |
256 |
|
|
SemiNaiveReduced(ires+(m*size), |
257 |
|
|
bw, |
258 |
|
|
m, |
259 |
|
|
fltres, |
260 |
|
|
scratchpad, |
261 |
|
|
seminaive_naive_table[m], |
262 |
|
|
weights, |
263 |
|
|
dctPlan); |
264 |
|
|
|
265 |
|
|
/* now load imaginary part of coefficients into output space */ |
266 |
|
|
memcpy(idataptr, fltres, sizeof(double) * (bw - m)); |
267 |
|
|
|
268 |
|
|
idataptr += bw-m; |
269 |
|
|
|
270 |
|
|
} |
271 |
|
|
else{ |
272 |
|
|
/* do real part */ |
273 |
|
|
Naive_AnalysisX(rres+(m*size), |
274 |
|
|
bw, |
275 |
|
|
m, |
276 |
|
|
weights, |
277 |
|
|
fltres, |
278 |
|
|
seminaive_naive_table[m], |
279 |
|
|
scratchpad ); |
280 |
|
|
memcpy(rdataptr, fltres, sizeof(double) * (bw - m)); |
281 |
|
|
rdataptr += bw-m; |
282 |
|
|
|
283 |
|
|
/* do imaginary part */ |
284 |
|
|
Naive_AnalysisX(ires+(m*size), |
285 |
|
|
bw, |
286 |
|
|
m, |
287 |
|
|
weights, |
288 |
|
|
fltres, |
289 |
|
|
seminaive_naive_table[m], |
290 |
|
|
scratchpad ); |
291 |
|
|
memcpy(idataptr, fltres, sizeof(double) * (bw - m)); |
292 |
|
|
idataptr += bw-m; |
293 |
|
|
} |
294 |
|
|
} |
295 |
|
|
|
296 |
|
|
/*** now do upper coefficients ****/ |
297 |
|
|
|
298 |
|
|
/* now if the data is real, we don't have to compute the |
299 |
|
|
coefficients whose order is less than 0, i.e. since |
300 |
|
|
the data is real, we know that |
301 |
|
|
f-hat(l,-m) = (-1)^m * conjugate(f-hat(l,m)), |
302 |
|
|
so use that to get the rest of the coefficients |
303 |
|
|
|
304 |
|
|
dataformat =0 -> samples are complex, =1 -> samples real |
305 |
|
|
|
306 |
|
|
*/ |
307 |
|
|
|
308 |
|
|
if( dataformat == 0 ){ |
309 |
|
|
|
310 |
|
|
/* note that m is greater than bw here, but this is for |
311 |
|
|
purposes of indexing the input data arrays. |
312 |
|
|
The "true" value of m as a parameter for Pml is |
313 |
|
|
size - m */ |
314 |
|
|
|
315 |
|
|
for (m=bw+1; m<size; m++) { |
316 |
|
|
/* |
317 |
|
|
fprintf(stderr,"m = %d\n",-(size-m)); |
318 |
|
|
*/ |
319 |
|
|
|
320 |
|
|
if ( (size-m) < cutoff ) |
321 |
|
|
{ |
322 |
|
|
/* do real part */ |
323 |
|
|
SemiNaiveReduced(rres+(m*size), |
324 |
|
|
bw, |
325 |
|
|
size-m, |
326 |
|
|
fltres, |
327 |
|
|
scratchpad, |
328 |
|
|
seminaive_naive_table[size-m], |
329 |
|
|
weights, |
330 |
|
|
dctPlan ); |
331 |
|
|
|
332 |
|
|
/* now load real part of coefficients into output space */ |
333 |
|
|
if ((m % 2) != 0) { |
334 |
|
|
for (i=0; i<m-bw; i++) |
335 |
|
|
rdataptr[i] = -fltres[i]; |
336 |
|
|
} |
337 |
|
|
else { |
338 |
|
|
memcpy(rdataptr, fltres, sizeof(double) * (m - bw)); |
339 |
|
|
} |
340 |
|
|
rdataptr += m-bw; |
341 |
|
|
|
342 |
|
|
/* do imaginary part */ |
343 |
|
|
SemiNaiveReduced(ires+(m*size), |
344 |
|
|
bw, |
345 |
|
|
size-m, |
346 |
|
|
fltres, |
347 |
|
|
scratchpad, |
348 |
|
|
seminaive_naive_table[size-m], |
349 |
|
|
weights, |
350 |
|
|
dctPlan); |
351 |
|
|
|
352 |
|
|
/* now load imag part of coefficients into output space */ |
353 |
|
|
if ((m % 2) != 0) { |
354 |
|
|
for (i=0; i<m-bw; i++) |
355 |
|
|
idataptr[i] = -fltres[i]; |
356 |
|
|
} |
357 |
|
|
else { |
358 |
|
|
memcpy(idataptr, fltres, sizeof(double) * (m - bw)); |
359 |
|
|
} |
360 |
|
|
idataptr += m-bw; |
361 |
|
|
} |
362 |
|
|
else |
363 |
|
|
{ |
364 |
|
|
Naive_AnalysisX(rres+(m*size), |
365 |
|
|
bw, |
366 |
|
|
size-m, |
367 |
|
|
weights, |
368 |
|
|
fltres, |
369 |
|
|
seminaive_naive_table[size-m], |
370 |
|
|
scratchpad); |
371 |
|
|
|
372 |
|
|
/* now load real part of coefficients into output space */ |
373 |
|
|
if ((m % 2) != 0) { |
374 |
|
|
for (i=0; i<m-bw; i++) |
375 |
|
|
rdataptr[i] = -fltres[i]; |
376 |
|
|
} |
377 |
|
|
else { |
378 |
|
|
memcpy(rdataptr, fltres, sizeof(double) * (m - bw)); |
379 |
|
|
} |
380 |
|
|
rdataptr += m-bw; |
381 |
|
|
|
382 |
|
|
/* do imaginary part */ |
383 |
|
|
Naive_AnalysisX(ires+(m*size), |
384 |
|
|
bw, |
385 |
|
|
size-m, |
386 |
|
|
weights, |
387 |
|
|
fltres, |
388 |
|
|
seminaive_naive_table[size-m], |
389 |
|
|
scratchpad); |
390 |
|
|
|
391 |
|
|
/* now load imag part of coefficients into output space */ |
392 |
|
|
if ((m % 2) != 0) { |
393 |
|
|
for (i=0; i<m-bw; i++) |
394 |
|
|
idataptr[i] = -fltres[i]; |
395 |
|
|
} |
396 |
|
|
else { |
397 |
|
|
memcpy(idataptr, fltres, sizeof(double) * (m - bw)); |
398 |
|
|
} |
399 |
|
|
idataptr += m-bw; |
400 |
|
|
|
401 |
|
|
} |
402 |
|
|
|
403 |
|
|
} |
404 |
|
|
} |
405 |
|
|
else /**** if the data is real ****/ |
406 |
|
|
{ |
407 |
|
|
pow_one = 1.0; |
408 |
|
|
for(i = 1; i < bw; i++){ |
409 |
|
|
pow_one *= -1.0; |
410 |
|
|
for( j = i; j < bw; j++){ |
411 |
|
|
/* |
412 |
|
|
SEANINDEXP(dummy0,i,j,bw); |
413 |
|
|
SEANINDEXN(dummy1,-i,j,bw); |
414 |
|
|
*/ |
415 |
|
|
|
416 |
|
|
dummy0 = seanindex(i, j, bw); |
417 |
|
|
dummy1 = seanindex(-i, j, bw); |
418 |
|
|
|
419 |
|
|
rcoeffs[dummy1] = |
420 |
|
|
pow_one * rcoeffs[dummy0]; |
421 |
|
|
icoeffs[dummy1] = |
422 |
|
|
-pow_one * icoeffs[dummy0]; |
423 |
|
|
} |
424 |
|
|
} |
425 |
|
|
} |
426 |
|
|
|
427 |
|
|
} |
428 |
|
|
|