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 |
|
|
/* source code for generating cosine transforms |
39 |
|
|
of Pml and Gml functions */ |
40 |
|
|
|
41 |
|
|
#include <math.h> |
42 |
|
|
#include <string.h> /* to declare memcpy */ |
43 |
|
|
#include <stdlib.h> |
44 |
|
|
#include <stdio.h> |
45 |
|
|
|
46 |
|
|
#include "fftw3.h" |
47 |
|
|
#include "primitive.h" |
48 |
|
|
#include "pmls.h" |
49 |
|
|
|
50 |
|
|
|
51 |
|
|
/************************************************************************/ |
52 |
|
|
/* utility functions for table management */ |
53 |
|
|
/************************************************************************/ |
54 |
|
|
/* Computes the number of non-zero entries in a table containing |
55 |
|
|
cosine series coefficients of all the P(m,l) or G(m,l) functions |
56 |
|
|
necessary for computing the seminaive transform for a given |
57 |
|
|
bandwidth bw and order m. Works specifically for tables |
58 |
|
|
generated by CosPmlTableGen() |
59 |
|
|
*/ |
60 |
|
|
|
61 |
|
|
int TableSize(int m, |
62 |
|
|
int bw) |
63 |
|
|
{ |
64 |
|
|
|
65 |
|
|
int k ; |
66 |
|
|
int fudge, fudge2 ; |
67 |
|
|
int a1, a2, a3 ; |
68 |
|
|
|
69 |
|
|
if ( bw % 2 ) /* if the bandwidth is odd */ |
70 |
|
|
{ |
71 |
|
|
k = bw/2 ; |
72 |
|
|
fudge = (m+1)%2 ; |
73 |
|
|
|
74 |
|
|
a1 = k*(k+1); |
75 |
|
|
a2 = fudge*(k+1); |
76 |
|
|
|
77 |
|
|
fudge2 = m/2; |
78 |
|
|
a3 = fudge2*(fudge2+1); |
79 |
|
|
|
80 |
|
|
} |
81 |
|
|
else /* bandwidth is even */ |
82 |
|
|
{ |
83 |
|
|
k = bw/2 ; |
84 |
|
|
fudge = m%2 ; |
85 |
|
|
|
86 |
|
|
a1 = (k-fudge)*(k-fudge+1); |
87 |
|
|
a2 = fudge*k; |
88 |
|
|
|
89 |
|
|
fudge2 = m/2; |
90 |
|
|
a3 = fudge2*(fudge2+1); |
91 |
|
|
|
92 |
|
|
} |
93 |
|
|
|
94 |
|
|
return(a1+a2-a3); |
95 |
|
|
|
96 |
|
|
} |
97 |
|
|
/************************************************************************/ |
98 |
|
|
/* Spharmonic_TableSize(bw) returns an integer value for |
99 |
|
|
the amount of space necessary to fill out an entire spharmonic |
100 |
|
|
table. Note that in the above TableSize() formula, |
101 |
|
|
you need to sum this formula over m as m ranges from 0 to |
102 |
|
|
(bw-1). The critical closed form that you need is that |
103 |
|
|
|
104 |
|
|
\sum_{k=0}^n = \frac{(n(n+1)(2n+1)}{6} |
105 |
|
|
|
106 |
|
|
You also need to account for integer division. |
107 |
|
|
From this you should derive an upper bound on the |
108 |
|
|
amount of space. |
109 |
|
|
|
110 |
|
|
Some notes - because of integer division, you need to account for |
111 |
|
|
a fudge factor - this is the additional " + bw" at the |
112 |
|
|
end. This gaurantees that you will always have slightly more |
113 |
|
|
space than you need, which is clearly better than underestimating! |
114 |
|
|
Also, if bw > 512, the closed form |
115 |
|
|
fails because of the bw*bw*bw term (at least on Sun Sparcstations) |
116 |
|
|
so the loop computation is used instead. |
117 |
|
|
|
118 |
|
|
Also, the transpose is exactly the same size, obviously. |
119 |
|
|
|
120 |
|
|
*/ |
121 |
|
|
|
122 |
|
|
int Spharmonic_TableSize(int bw) |
123 |
|
|
{ |
124 |
|
|
int m, sum; |
125 |
|
|
|
126 |
|
|
if (bw > 512) |
127 |
|
|
{ |
128 |
|
|
sum = 0; |
129 |
|
|
|
130 |
|
|
for (m=0; m<bw; m++) |
131 |
|
|
sum += TableSize(m,bw); |
132 |
|
|
|
133 |
|
|
return sum; |
134 |
|
|
} |
135 |
|
|
else |
136 |
|
|
{ |
137 |
|
|
return ( |
138 |
|
|
(((4*(bw*bw*bw)) + (6*(bw*bw)) - (8*bw))/24) |
139 |
|
|
+ bw |
140 |
|
|
); |
141 |
|
|
} |
142 |
|
|
} |
143 |
|
|
|
144 |
|
|
/************************************************************************/ |
145 |
|
|
/* Reduced_Spharmonic_TableSize(bw,m) returns an integer value for |
146 |
|
|
the amount of space necessary to fill out a spharmonic table |
147 |
|
|
if interesting in using it only for orders up to (but NOT |
148 |
|
|
including) order m. |
149 |
|
|
This will be used in the hybrid algorithm's call of the |
150 |
|
|
semi-naive algorithm (which won't need the full table ... hopefully |
151 |
|
|
this'll cut down on the memory usage). |
152 |
|
|
|
153 |
|
|
Also, the transpose is exactly the same size, obviously. |
154 |
|
|
|
155 |
|
|
This is a "reduced" version of Spharmonic_TableSize(m). |
156 |
|
|
|
157 |
|
|
*/ |
158 |
|
|
|
159 |
|
|
int Reduced_SpharmonicTableSize(int bw, |
160 |
|
|
int m) |
161 |
|
|
{ |
162 |
|
|
|
163 |
|
|
int i, sum; |
164 |
|
|
|
165 |
|
|
sum = 0; |
166 |
|
|
|
167 |
|
|
for (i=0; i<m; i++) |
168 |
|
|
sum += TableSize(i,bw); |
169 |
|
|
|
170 |
|
|
return sum; |
171 |
|
|
} |
172 |
|
|
|
173 |
|
|
|
174 |
|
|
/************************************************************************/ |
175 |
|
|
/* For an array containing cosine series coefficients of Pml or Gml |
176 |
|
|
functions, computes the location of the first coefficient of Pml. |
177 |
|
|
This supersedes the TableOffset() function. |
178 |
|
|
Assumes table is generated by CosPmlTableGen() |
179 |
|
|
*/ |
180 |
|
|
|
181 |
|
|
int NewTableOffset(int m, |
182 |
|
|
int l) |
183 |
|
|
{ |
184 |
|
|
int offset; |
185 |
|
|
int tm, tl; |
186 |
|
|
|
187 |
|
|
if ( m % 2 ) |
188 |
|
|
{ |
189 |
|
|
tl = l-1; |
190 |
|
|
tm = m-1; |
191 |
|
|
} |
192 |
|
|
else |
193 |
|
|
{ |
194 |
|
|
tl = l; |
195 |
|
|
tm = m; |
196 |
|
|
} |
197 |
|
|
|
198 |
|
|
offset = ((tl/2)*((tl/2)+1)) - ((tm/2)*((tm/2)+1)); |
199 |
|
|
if (tl % 2) |
200 |
|
|
offset += (tl/2)+1; |
201 |
|
|
|
202 |
|
|
return offset; |
203 |
|
|
} |
204 |
|
|
|
205 |
|
|
/************************************************************************/ |
206 |
|
|
/* |
207 |
|
|
generate all of the cosine series for L2-normalized Pmls or Gmls for |
208 |
|
|
a specified value of m. Note especially that since series are |
209 |
|
|
zero-striped, all zeroes have been removed. |
210 |
|
|
|
211 |
|
|
tablespace points to a double array of size TableSize(m,bw); |
212 |
|
|
|
213 |
|
|
Workspace needs to be |
214 |
|
|
9 * bw |
215 |
|
|
|
216 |
|
|
Let P(m,l,j) represent the jth coefficient of the |
217 |
|
|
cosine series representation of Pml. The array |
218 |
|
|
stuffed into tablespace is organized as follows: |
219 |
|
|
|
220 |
|
|
P(m,m,0) P(m,m,2) ... P(m,m,m) |
221 |
|
|
P(m,m+1,1) P(m,m+1,3)... P(m,m+1,m+1) |
222 |
|
|
P(m,m+2,0) P(m,m+2,2) ... P(m,m+2,m+2) |
223 |
|
|
|
224 |
|
|
etc. Appropriate modifications are made for m odd (Gml functions). |
225 |
|
|
|
226 |
|
|
|
227 |
|
|
NOTE that the Pmls or Gmls are being sampled at bw-many points, |
228 |
|
|
and not 2*bw-many points. I can get away with this. HOWEVER, I |
229 |
|
|
need to multiply the coefficients by sqrt(2), because the expected |
230 |
|
|
input of the seminaive transform of bandwidth bw will be sampled |
231 |
|
|
at 2-bw many points. So the sqrt(2) is a scaling factor. |
232 |
|
|
|
233 |
|
|
|
234 |
|
|
*/ |
235 |
|
|
|
236 |
|
|
void CosPmlTableGen(int bw, |
237 |
|
|
int m, |
238 |
|
|
double *tablespace, |
239 |
|
|
double *workspace) |
240 |
|
|
{ |
241 |
|
|
double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4; |
242 |
|
|
double *x_i, *eval_args; |
243 |
|
|
double *tableptr, *cosres ; |
244 |
|
|
int i, j, k; |
245 |
|
|
|
246 |
|
|
/* fftw stuff now */ |
247 |
|
|
double fudge ; |
248 |
|
|
fftw_plan p ; |
249 |
|
|
|
250 |
|
|
prevprev = workspace; |
251 |
|
|
prev = prevprev + bw; |
252 |
|
|
temp1 = prev + bw; |
253 |
|
|
temp2 = temp1 + bw; |
254 |
|
|
temp3 = temp2 + bw; |
255 |
|
|
temp4 = temp3 + bw; |
256 |
|
|
x_i = temp4 + bw; |
257 |
|
|
eval_args = x_i + bw; |
258 |
|
|
cosres = eval_args + bw; |
259 |
|
|
|
260 |
|
|
tableptr = tablespace; |
261 |
|
|
|
262 |
|
|
/* make fftw plan */ |
263 |
|
|
p = fftw_plan_r2r_1d( bw, temp4, cosres, |
264 |
|
|
FFTW_REDFT10, FFTW_ESTIMATE ) ; |
265 |
|
|
|
266 |
|
|
/* main loop */ |
267 |
|
|
|
268 |
|
|
/* Set the initial number of evaluation points to appropriate |
269 |
|
|
amount */ |
270 |
|
|
|
271 |
|
|
/* now get the evaluation nodes */ |
272 |
|
|
EvalPts(bw,x_i); |
273 |
|
|
ArcCosEvalPts(bw,eval_args); |
274 |
|
|
|
275 |
|
|
/* set initial values of first two Pmls */ |
276 |
|
|
for (i=0; i<bw; i++) |
277 |
|
|
prevprev[i] = 0.0; |
278 |
|
|
|
279 |
|
|
if (m == 0) |
280 |
|
|
for (i=0; i<bw; i++) |
281 |
|
|
prev[i] = 0.707106781186547; /* sqrt(1/2) */ |
282 |
|
|
else |
283 |
|
|
Pmm_L2(m, eval_args, bw, prev); |
284 |
|
|
|
285 |
|
|
|
286 |
|
|
if ( m % 2 ) /* need to divide out sin x */ |
287 |
|
|
for (i=0; i<bw; i++) |
288 |
|
|
prev[i] /= sin(eval_args[i]); |
289 |
|
|
|
290 |
|
|
|
291 |
|
|
/* set k to highest degree coefficient */ |
292 |
|
|
if ((m % 2) == 0) |
293 |
|
|
k = m; |
294 |
|
|
else |
295 |
|
|
k = m-1; |
296 |
|
|
|
297 |
|
|
/* now compute cosine transform */ |
298 |
|
|
memcpy( temp4, prev, sizeof(double) * bw ); |
299 |
|
|
fftw_execute( p ); |
300 |
|
|
cosres[0] *= 0.707106781186547 ; |
301 |
|
|
fudge = 1. / sqrt(((double) bw ) ); |
302 |
|
|
for ( i = 0 ; i < bw ; i ++ ) |
303 |
|
|
cosres[i] *= fudge ; |
304 |
|
|
|
305 |
|
|
/* store what I've got so far */ |
306 |
|
|
for (i=0; i<=k; i+=2) |
307 |
|
|
tableptr[i/2] = cosres[i]; |
308 |
|
|
|
309 |
|
|
/* update tableptr */ |
310 |
|
|
tableptr += k/2+1; |
311 |
|
|
|
312 |
|
|
/* now generate remaining pmls */ |
313 |
|
|
for (i=0; i<bw-m-1; i++) |
314 |
|
|
{ |
315 |
|
|
vec_mul(L2_cn(m,m+i),prevprev,temp1,bw); |
316 |
|
|
vec_pt_mul(prev, x_i, temp2, bw); |
317 |
|
|
vec_mul(L2_an(m,m+i), temp2, temp3, bw); |
318 |
|
|
vec_add(temp3, temp1, temp4, bw); /* temp4 now contains P(m,m+i+1) */ |
319 |
|
|
|
320 |
|
|
/* compute cosine transform */ |
321 |
|
|
fftw_execute( p ); |
322 |
|
|
cosres[0] *= 0.707106781186547 ; |
323 |
|
|
for ( j = 0 ; j < bw ; j ++ ) |
324 |
|
|
cosres[j] *= fudge ; |
325 |
|
|
|
326 |
|
|
/* update degree counter */ |
327 |
|
|
k++; |
328 |
|
|
|
329 |
|
|
/* now put decimated result into table */ |
330 |
|
|
if ( i % 2 ) |
331 |
|
|
for (j=0; j<=k; j+=2) |
332 |
|
|
tableptr[j/2] = cosres[j]; |
333 |
|
|
else |
334 |
|
|
for (j=1; j<=k; j+=2) |
335 |
|
|
tableptr[j/2] = cosres[j]; |
336 |
|
|
|
337 |
|
|
/* update tableptr */ |
338 |
|
|
tableptr += k/2+1; |
339 |
|
|
|
340 |
|
|
/* now update Pi and P(i+1) */ |
341 |
|
|
memcpy(prevprev, prev, sizeof(double) * bw); |
342 |
|
|
memcpy(prev, temp4, sizeof(double) * bw); |
343 |
|
|
} |
344 |
|
|
|
345 |
|
|
fftw_destroy_plan( p ); |
346 |
|
|
|
347 |
|
|
} |
348 |
|
|
|
349 |
|
|
|
350 |
|
|
/************************************************************************/ |
351 |
|
|
/* RowSize returns the number of non-zero coefficients in a row of the |
352 |
|
|
cospmltable if were really in matrix form. Helpful in transpose |
353 |
|
|
computations. It is helpful to think of the parameter l as |
354 |
|
|
the row of the corresponding matrix. |
355 |
|
|
*/ |
356 |
|
|
|
357 |
|
|
int RowSize(int m, |
358 |
|
|
int l) |
359 |
|
|
{ |
360 |
|
|
if (l < m) |
361 |
|
|
return 0; |
362 |
|
|
else |
363 |
|
|
{ |
364 |
|
|
if ((m % 2) == 0) |
365 |
|
|
return ((l/2)+1); |
366 |
|
|
else |
367 |
|
|
return (((l-1)/2)+1); |
368 |
|
|
} |
369 |
|
|
} |
370 |
|
|
/************************************************************************/ |
371 |
|
|
/* Transposed row size returns the number of non-zero coefficients |
372 |
|
|
in the transposition of the matrix representing a cospmltable. |
373 |
|
|
Used for generating arrays for inverse seminaive transform. |
374 |
|
|
Unlike RowSize, need to know the bandwidth bw. Also, in |
375 |
|
|
the cospml array, the first m+1 rows are empty, but in |
376 |
|
|
the transpose, all rows have non-zero entries, and the first |
377 |
|
|
m+1 columns are empty. So the input parameters are a bit different |
378 |
|
|
in the you need to specify the row you want. |
379 |
|
|
|
380 |
|
|
*/ |
381 |
|
|
|
382 |
|
|
int Transpose_RowSize(int row, |
383 |
|
|
int m, |
384 |
|
|
int bw) |
385 |
|
|
{ |
386 |
|
|
/* my version might be longer, but at least I understand |
387 |
|
|
it better, and it's only minimally recursive */ |
388 |
|
|
|
389 |
|
|
if ( bw % 2 ) |
390 |
|
|
{ |
391 |
|
|
if ( m % 2 ) |
392 |
|
|
{ |
393 |
|
|
if ( m == 1 ) |
394 |
|
|
return( (bw-row)/2 ); |
395 |
|
|
else if ( row < m - 1 ) |
396 |
|
|
return ( (bw-m+1)/2 ); |
397 |
|
|
else |
398 |
|
|
return ( Transpose_RowSize(row, 1, bw) ) ; |
399 |
|
|
} |
400 |
|
|
else |
401 |
|
|
{ |
402 |
|
|
if ( m == 0 ) |
403 |
|
|
return( (bw-row)/2 + ((row+1)%2) ); |
404 |
|
|
else if ( row < m ) |
405 |
|
|
return ( (bw-m)/2 + ((row+1)%2) ); |
406 |
|
|
else |
407 |
|
|
return ( Transpose_RowSize(row, 0, bw) ) ; |
408 |
|
|
} |
409 |
|
|
} |
410 |
|
|
else |
411 |
|
|
{ |
412 |
|
|
if ( m % 2 ) |
413 |
|
|
{ |
414 |
|
|
if ( m == 1 ) |
415 |
|
|
return( (bw-row)/2 ); |
416 |
|
|
else if ( row < m - 1 ) |
417 |
|
|
return ( (bw-m+1)/2 - (row%2) ); |
418 |
|
|
else |
419 |
|
|
return ( Transpose_RowSize(row, 1, bw) ) ; |
420 |
|
|
} |
421 |
|
|
else |
422 |
|
|
{ |
423 |
|
|
if ( m == 0 ) |
424 |
|
|
return( (bw-row)/2 + (row%2) ); |
425 |
|
|
else if ( row < m ) |
426 |
|
|
return ( (bw-m)/2 ); |
427 |
|
|
else |
428 |
|
|
return ( Transpose_RowSize(row, 0, bw) ) ; |
429 |
|
|
} |
430 |
|
|
} |
431 |
|
|
|
432 |
|
|
|
433 |
|
|
|
434 |
|
|
/*** original version |
435 |
|
|
|
436 |
|
|
if (row >= bw) |
437 |
|
|
return 0; |
438 |
|
|
else if ((m % 2) == 0) |
439 |
|
|
{ |
440 |
|
|
if (row <= m) |
441 |
|
|
return ( ((bw-m)/2) ); |
442 |
|
|
else |
443 |
|
|
return ( ((bw-row-1)/2) + 1); |
444 |
|
|
} |
445 |
|
|
else |
446 |
|
|
{ |
447 |
|
|
if (row == (bw-1)) |
448 |
|
|
return 0; |
449 |
|
|
else if (row >= m) |
450 |
|
|
return (Transpose_RowSize(row+1,m-1,bw)); |
451 |
|
|
else |
452 |
|
|
return (Transpose_RowSize(row+1,m-1,bw) - (row % 2)); |
453 |
|
|
} |
454 |
|
|
|
455 |
|
|
***/ |
456 |
|
|
|
457 |
|
|
} |
458 |
|
|
|
459 |
|
|
/************************************************************************/ |
460 |
|
|
/* Inverse transform is transposition of forward transform. |
461 |
|
|
Thus, need to provide transposed version of table |
462 |
|
|
returned by CosPmlTableGen. This function does that |
463 |
|
|
by taking as input a cos_pml_table for a particular value |
464 |
|
|
of bw and m, and loads the result as a |
465 |
|
|
transposed, decimated version of it for use by an inverse |
466 |
|
|
seminaive transform computation. |
467 |
|
|
|
468 |
|
|
result needs to be of size TableSize(m,bw) |
469 |
|
|
|
470 |
|
|
*/ |
471 |
|
|
|
472 |
|
|
void Transpose_CosPmlTableGen(int bw, |
473 |
|
|
int m, |
474 |
|
|
double *cos_pml_table, |
475 |
|
|
double *result) |
476 |
|
|
{ |
477 |
|
|
/* recall that cospml_table has had all the zeroes |
478 |
|
|
stripped out, and that if m is odd, then it is |
479 |
|
|
really a Gml function, which affects indexing a bit. |
480 |
|
|
*/ |
481 |
|
|
|
482 |
|
|
double *trans_tableptr, *tableptr; |
483 |
|
|
int i, row, rowsize, stride, offset, costable_offset; |
484 |
|
|
|
485 |
|
|
/* note that the number of non-zero entries is the same |
486 |
|
|
as in the non-transposed case */ |
487 |
|
|
|
488 |
|
|
trans_tableptr = result; |
489 |
|
|
|
490 |
|
|
/* now traverse the cos_pml_table , loading appropriate values |
491 |
|
|
into the rows of transposed array */ |
492 |
|
|
|
493 |
|
|
if ( m == bw - 1 ) |
494 |
|
|
memcpy( result, cos_pml_table, sizeof(double)*TableSize(m,bw)); |
495 |
|
|
else |
496 |
|
|
{ |
497 |
|
|
|
498 |
|
|
for (row = 0; row < bw; row++) |
499 |
|
|
{ |
500 |
|
|
/* if m odd, no need to do last row - all zeroes */ |
501 |
|
|
if (row == (bw-1)) |
502 |
|
|
{ |
503 |
|
|
if ( m % 2 ) |
504 |
|
|
return; |
505 |
|
|
} |
506 |
|
|
|
507 |
|
|
/* get the rowsize for the transposed array */ |
508 |
|
|
rowsize = Transpose_RowSize(row, m, bw); |
509 |
|
|
|
510 |
|
|
/* compute the starting point for values in cos_pml_table */ |
511 |
|
|
if (row <= m) |
512 |
|
|
{ |
513 |
|
|
if ((row % 2) == 0) |
514 |
|
|
tableptr = cos_pml_table + (row/2); |
515 |
|
|
else |
516 |
|
|
tableptr = cos_pml_table + (m/2) + 1 + (row/2); |
517 |
|
|
} |
518 |
|
|
else |
519 |
|
|
{ |
520 |
|
|
/* if row > m, then the highest degree coefficient |
521 |
|
|
of P(m,row) should be the first coefficient loaded |
522 |
|
|
into the transposed array, so figure out where |
523 |
|
|
this point is. |
524 |
|
|
*/ |
525 |
|
|
offset = 0; |
526 |
|
|
if ( (m%2) == 0 ) |
527 |
|
|
{ |
528 |
|
|
for (i=m; i<=row; i++) |
529 |
|
|
offset += RowSize(m, i); |
530 |
|
|
} |
531 |
|
|
else |
532 |
|
|
{ |
533 |
|
|
for (i=m;i<=row+1;i++) |
534 |
|
|
offset += RowSize(m, i); |
535 |
|
|
} |
536 |
|
|
/* now we are pointing one element too far, so decrement */ |
537 |
|
|
offset--; |
538 |
|
|
|
539 |
|
|
tableptr = cos_pml_table + offset; |
540 |
|
|
} |
541 |
|
|
|
542 |
|
|
/* stride is how far we need to jump between |
543 |
|
|
values in cos_pml_table, i.e., to traverse the columns of the |
544 |
|
|
cos_pml_table. Need to set initial value. Stride always |
545 |
|
|
increases by 2 after that |
546 |
|
|
*/ |
547 |
|
|
if (row <= m) |
548 |
|
|
stride = m + 2 - (m % 2) + (row % 2); |
549 |
|
|
else |
550 |
|
|
stride = row + 2; |
551 |
|
|
|
552 |
|
|
/* now load up this row of the transposed table */ |
553 |
|
|
costable_offset = 0; |
554 |
|
|
for (i=0; i < rowsize; i++) |
555 |
|
|
{ |
556 |
|
|
trans_tableptr[i] = tableptr[costable_offset]; |
557 |
|
|
costable_offset += stride; |
558 |
|
|
stride += 2; |
559 |
|
|
|
560 |
|
|
} /* closes i loop */ |
561 |
|
|
|
562 |
|
|
trans_tableptr += rowsize; |
563 |
|
|
|
564 |
|
|
} /* closes row loop */ |
565 |
|
|
} |
566 |
|
|
|
567 |
|
|
} |
568 |
|
|
/************************************************************************/ |
569 |
|
|
/* This is a function that returns all of the (cosine transforms of) |
570 |
|
|
Pmls and Gmls necessary |
571 |
|
|
to do a full spherical harmonic transform, i.e., it calls |
572 |
|
|
CosPmlTableGen for each value of m less than bw, returning a |
573 |
|
|
table of tables ( a pointer of type (double **), which points |
574 |
|
|
to an array of size m, each containing a (double *) pointer |
575 |
|
|
to a set of cospml or cosgml values, which are the (decimated) |
576 |
|
|
cosine series representations of Pml (even m) or Gml (odd m) |
577 |
|
|
functions. See CosPmlTableGen for further clarification. |
578 |
|
|
|
579 |
|
|
Inputs - the bandwidth bw of the problem |
580 |
|
|
resultspace - need to allocate Spharmonic_TableSize(bw) for storing results |
581 |
|
|
workspace - needs to be (16 * bw) |
582 |
|
|
|
583 |
|
|
Note that resultspace is necessary and contains the results/values |
584 |
|
|
so one should be careful about when it is OK to re-use this space. |
585 |
|
|
workspace, though, does not have any meaning after this function is |
586 |
|
|
finished executing |
587 |
|
|
|
588 |
|
|
*/ |
589 |
|
|
|
590 |
|
|
double **Spharmonic_Pml_Table(int bw, |
591 |
|
|
double *resultspace, |
592 |
|
|
double *workspace) |
593 |
|
|
{ |
594 |
|
|
|
595 |
|
|
int i; |
596 |
|
|
double **spharmonic_pml_table; |
597 |
|
|
|
598 |
|
|
/* allocate an array of double pointers */ |
599 |
|
|
spharmonic_pml_table = (double **) malloc(sizeof(double *) * bw); |
600 |
|
|
|
601 |
|
|
/* traverse the array, assigning a location in the resultspace |
602 |
|
|
to each pointer */ |
603 |
|
|
|
604 |
|
|
spharmonic_pml_table[0] = resultspace; |
605 |
|
|
|
606 |
|
|
for (i=1; i<bw; i++) |
607 |
|
|
{ |
608 |
|
|
spharmonic_pml_table[i] = spharmonic_pml_table[i-1] + |
609 |
|
|
TableSize(i-1,bw); |
610 |
|
|
} |
611 |
|
|
|
612 |
|
|
/* now load up the array with CosPml and CosGml values */ |
613 |
|
|
for (i=0; i<bw; i++) |
614 |
|
|
{ |
615 |
|
|
CosPmlTableGen(bw, i, spharmonic_pml_table[i], workspace); |
616 |
|
|
} |
617 |
|
|
|
618 |
|
|
/* that's it */ |
619 |
|
|
|
620 |
|
|
return spharmonic_pml_table; |
621 |
|
|
} |
622 |
|
|
|
623 |
|
|
|
624 |
|
|
/************************************************************************/ |
625 |
|
|
/* For the inverse semi-naive spharmonic transform, need the "transpose" |
626 |
|
|
of the spharmonic_pml_table. Need to be careful because the |
627 |
|
|
entries in the spharmonic_pml_table have been decimated, i.e., |
628 |
|
|
the zeroes have been stripped out. |
629 |
|
|
|
630 |
|
|
Inputs are a spharmonic_pml_table generated by Spharmonic_Pml_Table |
631 |
|
|
and the bandwidth bw |
632 |
|
|
|
633 |
|
|
Allocates memory for the (double **) result |
634 |
|
|
also allocates memory |
635 |
|
|
|
636 |
|
|
resultspace - need to allocate Spharmonic_TableSize(bw) for storing results |
637 |
|
|
workspace - not needed, but argument added to avoid |
638 |
|
|
confusion wth Spharmonic_Pml_Table |
639 |
|
|
|
640 |
|
|
*/ |
641 |
|
|
|
642 |
|
|
double **Transpose_Spharmonic_Pml_Table(double **spharmonic_pml_table, |
643 |
|
|
int bw, |
644 |
|
|
double *resultspace, |
645 |
|
|
double *workspace) |
646 |
|
|
{ |
647 |
|
|
|
648 |
|
|
int i; |
649 |
|
|
double **transpose_spharmonic_pml_table; |
650 |
|
|
|
651 |
|
|
/* allocate an array of double pointers */ |
652 |
|
|
transpose_spharmonic_pml_table = (double **) malloc(sizeof(double *) * bw); |
653 |
|
|
|
654 |
|
|
/* now need to load up the transpose_spharmonic_pml_table by transposing |
655 |
|
|
the tables in the spharmonic_pml_table */ |
656 |
|
|
|
657 |
|
|
transpose_spharmonic_pml_table[0] = resultspace; |
658 |
|
|
|
659 |
|
|
for (i=0; i<bw; i++) |
660 |
|
|
{ |
661 |
|
|
Transpose_CosPmlTableGen(bw, |
662 |
|
|
i, |
663 |
|
|
spharmonic_pml_table[i], |
664 |
|
|
transpose_spharmonic_pml_table[i]); |
665 |
|
|
|
666 |
|
|
if (i != (bw-1)) |
667 |
|
|
{ |
668 |
|
|
transpose_spharmonic_pml_table[i+1] = |
669 |
|
|
transpose_spharmonic_pml_table[i] + TableSize(i, bw); |
670 |
|
|
} |
671 |
|
|
} |
672 |
|
|
|
673 |
|
|
return transpose_spharmonic_pml_table; |
674 |
|
|
} |
675 |
|
|
|
676 |
|
|
|
677 |
|
|
/************************************************************************/ |
678 |
|
|
/* Reduced_Naive_TableSize(bw,m) returns an integer value for |
679 |
|
|
the amount of space necessary to fill out a reduced naive table |
680 |
|
|
of pmls if interested in using it only for orders m through bw - 1. |
681 |
|
|
|
682 |
|
|
*/ |
683 |
|
|
|
684 |
|
|
int Reduced_Naive_TableSize(int bw, |
685 |
|
|
int m) |
686 |
|
|
{ |
687 |
|
|
|
688 |
|
|
int i, sum; |
689 |
|
|
|
690 |
|
|
sum = 0; |
691 |
|
|
|
692 |
|
|
for (i=m; i<bw; i++) |
693 |
|
|
sum += ( 2 * bw * (bw - i)); |
694 |
|
|
|
695 |
|
|
return sum; |
696 |
|
|
|
697 |
|
|
} |
698 |
|
|
|
699 |
|
|
/************************************************************* |
700 |
|
|
|
701 |
|
|
just like Spharmonic_Pml_Table(), except generates a |
702 |
|
|
table for use with the semi-naive and naive algorithms. |
703 |
|
|
|
704 |
|
|
m is the cutoff order, where to switch from semi-naive to |
705 |
|
|
naive algorithms |
706 |
|
|
|
707 |
|
|
bw = bandwidth of problem |
708 |
|
|
m = where to switch algorithms (order where naive is FIRST done) |
709 |
|
|
resultspace = where to store results, must be of |
710 |
|
|
size Reduced_Naive_TableSize(bw, m) + |
711 |
|
|
Reduced_SpharmonicTableSize(bw, m); |
712 |
|
|
|
713 |
|
|
***********************************************************/ |
714 |
|
|
|
715 |
|
|
double **SemiNaive_Naive_Pml_Table(int bw, |
716 |
|
|
int m, |
717 |
|
|
double *resultspace, |
718 |
|
|
double *workspace) |
719 |
|
|
{ |
720 |
|
|
int i; |
721 |
|
|
double **seminaive_naive_table; |
722 |
|
|
int lastspace; |
723 |
|
|
|
724 |
|
|
seminaive_naive_table = (double **) malloc(sizeof(double) * (bw+1)); |
725 |
|
|
|
726 |
|
|
seminaive_naive_table[0] = resultspace; |
727 |
|
|
|
728 |
|
|
|
729 |
|
|
for (i=1; i<m; i++) |
730 |
|
|
{ |
731 |
|
|
seminaive_naive_table[i] = seminaive_naive_table[i - 1] + |
732 |
|
|
TableSize(i-1,bw); |
733 |
|
|
} |
734 |
|
|
|
735 |
|
|
if( m == 0) |
736 |
|
|
{ |
737 |
|
|
lastspace = 0; |
738 |
|
|
for (i=m+1; i<bw; i++) |
739 |
|
|
{ |
740 |
|
|
seminaive_naive_table[i] = seminaive_naive_table[i - 1] + |
741 |
|
|
(2 * bw * (bw - (i - 1))); |
742 |
|
|
} |
743 |
|
|
} |
744 |
|
|
else |
745 |
|
|
{ |
746 |
|
|
lastspace = TableSize(m-1,bw); |
747 |
|
|
seminaive_naive_table[m] = seminaive_naive_table[m-1] + |
748 |
|
|
lastspace; |
749 |
|
|
for (i=m+1; i<bw; i++) |
750 |
|
|
{ |
751 |
|
|
seminaive_naive_table[i] = seminaive_naive_table[i - 1] + |
752 |
|
|
(2 * bw * (bw - (i - 1))); |
753 |
|
|
} |
754 |
|
|
} |
755 |
|
|
|
756 |
|
|
/* now load up the array with CosPml and CosGml values */ |
757 |
|
|
for (i=0; i<m; i++) |
758 |
|
|
{ |
759 |
|
|
CosPmlTableGen(bw, i, seminaive_naive_table[i], workspace); |
760 |
|
|
} |
761 |
|
|
|
762 |
|
|
/* now load up pml values */ |
763 |
|
|
for(i=m; i<bw; i++) |
764 |
|
|
{ |
765 |
|
|
PmlTableGen(bw, i, seminaive_naive_table[i], workspace); |
766 |
|
|
} |
767 |
|
|
|
768 |
|
|
/* that's it */ |
769 |
|
|
|
770 |
|
|
return seminaive_naive_table; |
771 |
|
|
|
772 |
|
|
} |
773 |
|
|
|
774 |
|
|
|
775 |
|
|
/************************************************************************/ |
776 |
|
|
/* For the inverse seminaive_naive transform, need the "transpose" |
777 |
|
|
of the seminaive_naive_pml_table. Need to be careful because the |
778 |
|
|
entries in the seminaive portion have been decimated, i.e., |
779 |
|
|
the zeroes have been stripped out. |
780 |
|
|
|
781 |
|
|
Inputs are a seminaive_naive_pml_table generated by SemiNaive_Naive_Pml_Table |
782 |
|
|
and the bandwidth bw and cutoff order m |
783 |
|
|
|
784 |
|
|
Allocates memory for the (double **) result |
785 |
|
|
also allocates memory |
786 |
|
|
|
787 |
|
|
resultspace - need to allocate Reduced_Naive_TableSize(bw, m) + |
788 |
|
|
Reduced_SpharmonicTableSize(bw, m) for storing results |
789 |
|
|
workspace - size 16 * bw |
790 |
|
|
|
791 |
|
|
*/ |
792 |
|
|
|
793 |
|
|
double **Transpose_SemiNaive_Naive_Pml_Table(double **seminaive_naive_pml_table, |
794 |
|
|
int bw, |
795 |
|
|
int m, |
796 |
|
|
double *resultspace, |
797 |
|
|
double *workspace) |
798 |
|
|
{ |
799 |
|
|
|
800 |
|
|
int i, lastspace; |
801 |
|
|
double **trans_seminaive_naive_pml_table; |
802 |
|
|
|
803 |
|
|
/* allocate an array of double pointers */ |
804 |
|
|
trans_seminaive_naive_pml_table = (double **) malloc(sizeof(double *) * (bw+1)); |
805 |
|
|
|
806 |
|
|
/* now need to load up the transpose_seminaive_naive_pml_table by transposing |
807 |
|
|
the tables in the seminiave portion of seminaive_naive_pml_table */ |
808 |
|
|
|
809 |
|
|
trans_seminaive_naive_pml_table[0] = resultspace; |
810 |
|
|
|
811 |
|
|
|
812 |
|
|
for (i=1; i<m; i++) |
813 |
|
|
{ |
814 |
|
|
trans_seminaive_naive_pml_table[i] = |
815 |
|
|
trans_seminaive_naive_pml_table[i - 1] + |
816 |
|
|
TableSize(i-1,bw); |
817 |
|
|
} |
818 |
|
|
|
819 |
|
|
if( m == 0 ) |
820 |
|
|
{ |
821 |
|
|
lastspace = 0; |
822 |
|
|
for (i=m+1; i<bw; i++) |
823 |
|
|
{ |
824 |
|
|
trans_seminaive_naive_pml_table[i] = |
825 |
|
|
trans_seminaive_naive_pml_table[i - 1] + |
826 |
|
|
(2 * bw * (bw - (i - 1))); |
827 |
|
|
} |
828 |
|
|
} |
829 |
|
|
else |
830 |
|
|
{ |
831 |
|
|
lastspace = TableSize(m-1,bw); |
832 |
|
|
trans_seminaive_naive_pml_table[m] = |
833 |
|
|
trans_seminaive_naive_pml_table[m-1] + |
834 |
|
|
lastspace; |
835 |
|
|
|
836 |
|
|
for (i=m+1; i<bw; i++) |
837 |
|
|
{ |
838 |
|
|
trans_seminaive_naive_pml_table[i] = |
839 |
|
|
trans_seminaive_naive_pml_table[i - 1] + |
840 |
|
|
(2 * bw * (bw - (i - 1))); |
841 |
|
|
} |
842 |
|
|
} |
843 |
|
|
|
844 |
|
|
for (i=0; i<m; i++) |
845 |
|
|
{ |
846 |
|
|
Transpose_CosPmlTableGen(bw, |
847 |
|
|
i, |
848 |
|
|
seminaive_naive_pml_table[i], |
849 |
|
|
trans_seminaive_naive_pml_table[i]); |
850 |
|
|
|
851 |
|
|
if (i != (bw-1)) |
852 |
|
|
{ |
853 |
|
|
trans_seminaive_naive_pml_table[i+1] = |
854 |
|
|
trans_seminaive_naive_pml_table[i] + TableSize(i, bw); |
855 |
|
|
} |
856 |
|
|
} |
857 |
|
|
|
858 |
|
|
/* now load up pml values */ |
859 |
|
|
for(i=m; i<bw; i++) |
860 |
|
|
{ |
861 |
|
|
PmlTableGen(bw, i, trans_seminaive_naive_pml_table[i], workspace); |
862 |
|
|
} |
863 |
|
|
|
864 |
|
|
return trans_seminaive_naive_pml_table; |
865 |
|
|
} |