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 |
/******************************************************************** |
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 |
|