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 computing the Legendre transform where |
39 |
|
|
projections are carried out in cosine space, i.e., the |
40 |
|
|
"seminaive" algorithm. |
41 |
|
|
|
42 |
|
|
For a description, see the related paper or Sean's thesis. |
43 |
|
|
|
44 |
|
|
*/ |
45 |
|
|
|
46 |
|
|
#include <math.h> |
47 |
|
|
#include <stdio.h> |
48 |
|
|
#include <string.h> /** for memcpy **/ |
49 |
|
|
#include "fftw3.h" |
50 |
|
|
|
51 |
|
|
#include "cospmls.h" |
52 |
|
|
|
53 |
|
|
|
54 |
|
|
/************************************************************************/ |
55 |
|
|
/* InvSemiNaiveReduced computes the inverse Legendre transform |
56 |
|
|
using the transposed seminaive algorithm. Note that because |
57 |
|
|
the Legendre transform is orthogonal, the inverse can be |
58 |
|
|
computed by transposing the matrix formulation of the |
59 |
|
|
problem. |
60 |
|
|
|
61 |
|
|
The forward transform looks like |
62 |
|
|
|
63 |
|
|
l = PCWf |
64 |
|
|
|
65 |
|
|
where f is the data vector, W is a quadrature matrix, |
66 |
|
|
C is a cosine transform matrix, P is a matrix |
67 |
|
|
full of coefficients of the cosine series representation |
68 |
|
|
of each Pml function P(m,m) P(m,m+1) ... P(m,bw-1), |
69 |
|
|
and l is the (associated) Legendre series representation |
70 |
|
|
of f. |
71 |
|
|
|
72 |
|
|
So to do the inverse, you do |
73 |
|
|
|
74 |
|
|
f = trans(C) trans(P) l |
75 |
|
|
|
76 |
|
|
so you need to transpose the matrix P from the forward transform |
77 |
|
|
and then do a cosine series evaluation. No quadrature matrix |
78 |
|
|
is necessary. If order m is odd, then there is also a sin |
79 |
|
|
factor that needs to be accounted for. |
80 |
|
|
|
81 |
|
|
Note that this function was written to be part of a full |
82 |
|
|
spherical harmonic transform, so a lot of precomputation |
83 |
|
|
has been assumed. |
84 |
|
|
|
85 |
|
|
Input argument description |
86 |
|
|
|
87 |
|
|
coeffs - a double pointer to an array of length |
88 |
|
|
(bw-m) containing associated |
89 |
|
|
Legendre series coefficients. Assumed |
90 |
|
|
that first entry contains the P(m,m) |
91 |
|
|
coefficient. |
92 |
|
|
|
93 |
|
|
bw - problem bandwidth |
94 |
|
|
|
95 |
|
|
m - order of the associated Legendre functions |
96 |
|
|
|
97 |
|
|
result - a double pointer to an array of (2*bw) samples |
98 |
|
|
representing the evaluation of the Legendre |
99 |
|
|
series at (2*bw) Chebyshev nodes. |
100 |
|
|
|
101 |
|
|
trans_cos_pml_table - double pointer to array representing |
102 |
|
|
the linearized form of trans(P) above. |
103 |
|
|
See cospmls.{h,c} for a description |
104 |
|
|
of the function Transpose_CosPmlTableGen() |
105 |
|
|
which generates this array. |
106 |
|
|
|
107 |
|
|
sin_values - when m is odd, need to factor in the sin(x) that |
108 |
|
|
is factored out of the generation of the values |
109 |
|
|
in trans(P). |
110 |
|
|
|
111 |
|
|
workspace - a double array of size 2*bw -> temp space involving |
112 |
|
|
intermediate array |
113 |
|
|
|
114 |
|
|
fplan - pointer to fftw plan with input array being fcos |
115 |
|
|
and output being result; I'll probably be using the |
116 |
|
|
guru interface to execute - that way I can apply the |
117 |
|
|
same plan to different arrays; the plan should be |
118 |
|
|
|
119 |
|
|
fftw_plan_r2r_1d( 2*bw, fcos, result, |
120 |
|
|
FFTW_REDFT01, FFTW_ESTIMATE ); |
121 |
|
|
|
122 |
|
|
*/ |
123 |
|
|
|
124 |
|
|
void InvSemiNaiveReduced(double *coeffs, |
125 |
|
|
int bw, |
126 |
|
|
int m, |
127 |
|
|
double *result, |
128 |
|
|
double *trans_cos_pml_table, |
129 |
|
|
double *sin_values, |
130 |
|
|
double *workspace, |
131 |
|
|
fftw_plan *fplan ) |
132 |
|
|
{ |
133 |
|
|
double *trans_tableptr; |
134 |
|
|
double *assoc_offset; |
135 |
|
|
int i, j, rowsize; |
136 |
|
|
double *p; |
137 |
|
|
double *fcos, fcos0, fcos1, fcos2, fcos3; |
138 |
|
|
double fudge ; |
139 |
|
|
|
140 |
|
|
fcos = workspace ; |
141 |
|
|
|
142 |
|
|
/* for paranoia, zero out arrays */ |
143 |
|
|
memset( fcos, 0, sizeof(double) * 2 * bw ); |
144 |
|
|
memset( result, 0, sizeof(double) * 2 * bw ); |
145 |
|
|
|
146 |
|
|
trans_tableptr = trans_cos_pml_table; |
147 |
|
|
p = trans_cos_pml_table; |
148 |
|
|
|
149 |
|
|
/* main loop - compute each value of fcos |
150 |
|
|
|
151 |
|
|
Note that all zeroes have been stripped out of the |
152 |
|
|
trans_cos_pml_table, so indexing is somewhat complicated. |
153 |
|
|
*/ |
154 |
|
|
|
155 |
|
|
for (i=0; i<bw; i++) |
156 |
|
|
{ |
157 |
|
|
if (i == (bw-1)) |
158 |
|
|
{ |
159 |
|
|
if ( m % 2 ) |
160 |
|
|
{ |
161 |
|
|
fcos[bw-1] = 0.0; |
162 |
|
|
break; |
163 |
|
|
} |
164 |
|
|
} |
165 |
|
|
|
166 |
|
|
rowsize = Transpose_RowSize(i, m, bw); |
167 |
|
|
if (i > m) |
168 |
|
|
assoc_offset = coeffs + (i - m) + (m % 2); |
169 |
|
|
else |
170 |
|
|
assoc_offset = coeffs + (i % 2); |
171 |
|
|
|
172 |
|
|
fcos0 = 0.0 ; fcos1 = 0.0; fcos2 = 0.0; fcos3 = 0.0; |
173 |
|
|
|
174 |
|
|
for (j = 0; j < rowsize % 4; ++j) |
175 |
|
|
fcos0 += assoc_offset[2*j] * trans_tableptr[j]; |
176 |
|
|
|
177 |
|
|
for ( ; j < rowsize; j += 4){ |
178 |
|
|
fcos0 += assoc_offset[2*j] * trans_tableptr[j]; |
179 |
|
|
fcos1 += assoc_offset[2*(j+1)] * trans_tableptr[j+1]; |
180 |
|
|
fcos2 += assoc_offset[2*(j+2)] * trans_tableptr[j+2]; |
181 |
|
|
fcos3 += assoc_offset[2*(j+3)] * trans_tableptr[j+3]; |
182 |
|
|
} |
183 |
|
|
fcos[i] = fcos0 + fcos1 + fcos2 + fcos3 ; |
184 |
|
|
|
185 |
|
|
trans_tableptr += rowsize; |
186 |
|
|
} |
187 |
|
|
|
188 |
|
|
|
189 |
|
|
/* |
190 |
|
|
now we have the cosine series for the result, |
191 |
|
|
so now evaluate the cosine series at 2*bw Chebyshev nodes |
192 |
|
|
*/ |
193 |
|
|
|
194 |
|
|
/* scale coefficients prior to taking inverse DCT */ |
195 |
|
|
fudge = 0.5 / sqrt((double) bw) ; |
196 |
|
|
for ( j = 1 ; j < 2*bw ; j ++ ) |
197 |
|
|
fcos[j] *= fudge ; |
198 |
|
|
fcos[0] /= sqrt(2. * ((double) bw)); |
199 |
|
|
|
200 |
|
|
/* now take the inverse dct */ |
201 |
|
|
/* NOTE that I am using the guru interface */ |
202 |
|
|
fftw_execute_r2r( *fplan, |
203 |
|
|
fcos, result ); |
204 |
|
|
|
205 |
|
|
/* if m is odd, then need to multiply by sin(x) at Chebyshev nodes */ |
206 |
|
|
if ( m % 2 ) |
207 |
|
|
{ |
208 |
|
|
for (j=0; j<(2*bw); j++) |
209 |
|
|
result[j] *= sin_values[j]; |
210 |
|
|
} |
211 |
|
|
|
212 |
|
|
trans_tableptr = p; |
213 |
|
|
|
214 |
|
|
/* amscray */ |
215 |
|
|
|
216 |
|
|
} |
217 |
|
|
|
218 |
|
|
/************************************************************************/ |
219 |
|
|
|
220 |
|
|
/* SemiNaiveReduced computes the Legendre transform of data. |
221 |
|
|
This function has been designed to be a component in |
222 |
|
|
a full spherical harmonic transform. |
223 |
|
|
|
224 |
|
|
data - pointer to double array of size (2*bw) containing |
225 |
|
|
function to be transformed. Assumes sampling at Chebyshev nodes |
226 |
|
|
|
227 |
|
|
bw - bandwidth of the problem |
228 |
|
|
m - order of the problem. 0 <= m < bw |
229 |
|
|
|
230 |
|
|
result - pointer to double array of length bw for returning computed |
231 |
|
|
Legendre coefficients. Contains |
232 |
|
|
bw-m coeffs, with the <f,P(m,m)> coefficient |
233 |
|
|
located in result[0] |
234 |
|
|
|
235 |
|
|
cos_pml_table - a pointer to an array containing the cosine |
236 |
|
|
series coefficients of the Pmls (or Gmls) |
237 |
|
|
for this problem. This table can be computed |
238 |
|
|
using the CosPmlTableGen() function, and |
239 |
|
|
the offset for a particular Pml can be had |
240 |
|
|
by calling the function NewTableOffset(). |
241 |
|
|
The size of the table is computed using |
242 |
|
|
the TableSize() function. Note that |
243 |
|
|
since the cosine series are always zero-striped, |
244 |
|
|
the zeroes have been removed. |
245 |
|
|
|
246 |
|
|
weights -> ptr to double array of size 4*bw - this array holds |
247 |
|
|
the weights for both even (starting at weights[0]) |
248 |
|
|
and odd (weights[2*bw]) transforms |
249 |
|
|
|
250 |
|
|
|
251 |
|
|
workspace -> tmp space: ptr to double array of size 4*bw |
252 |
|
|
|
253 |
|
|
fplan -> pointer to fftw plan with input array being weighted_data |
254 |
|
|
and output being cos_data; I'll probably be using the |
255 |
|
|
guru interface to execute; the plan should be |
256 |
|
|
|
257 |
|
|
fftw_plan_r2r_1d( 2*bw, weighted_data, cos_data, |
258 |
|
|
FFTW_REDFT10, FFTW_ESTIMATE ) ; |
259 |
|
|
|
260 |
|
|
|
261 |
|
|
*/ |
262 |
|
|
|
263 |
|
|
void SemiNaiveReduced(double *data, |
264 |
|
|
int bw, |
265 |
|
|
int m, |
266 |
|
|
double *result, |
267 |
|
|
double *workspace, |
268 |
|
|
double *cos_pml_table, |
269 |
|
|
double *weights, |
270 |
|
|
fftw_plan *fplan ) |
271 |
|
|
{ |
272 |
|
|
int i, j, n; |
273 |
|
|
double result0, result1, result2, result3; |
274 |
|
|
double fudge ; |
275 |
|
|
double d_bw; |
276 |
|
|
int toggle ; |
277 |
|
|
double *pml_ptr, *weighted_data, *cos_data ; |
278 |
|
|
|
279 |
|
|
n = 2*bw; |
280 |
|
|
d_bw = (double) bw; |
281 |
|
|
|
282 |
|
|
weighted_data = workspace ; |
283 |
|
|
cos_data = weighted_data + (2*bw) ; |
284 |
|
|
|
285 |
|
|
/* for paranoia, zero out the result array */ |
286 |
|
|
memset( result, 0, sizeof(double)*(bw-m)); |
287 |
|
|
|
288 |
|
|
/* |
289 |
|
|
need to apply quadrature weights to the data and compute |
290 |
|
|
the cosine transform |
291 |
|
|
*/ |
292 |
|
|
if ( m % 2 ) |
293 |
|
|
for ( i = 0; i < n ; ++i ) |
294 |
|
|
weighted_data[i] = data[ i ] * weights[ 2*bw + i ]; |
295 |
|
|
else |
296 |
|
|
for ( i = 0; i < n ; ++i ) |
297 |
|
|
weighted_data[i] = data[ i ] * weights[ i ]; |
298 |
|
|
|
299 |
|
|
/* |
300 |
|
|
smooth the weighted signal |
301 |
|
|
*/ |
302 |
|
|
|
303 |
|
|
fftw_execute_r2r( *fplan, |
304 |
|
|
weighted_data, |
305 |
|
|
cos_data ); |
306 |
|
|
|
307 |
|
|
/* need to normalize */ |
308 |
|
|
cos_data[0] *= 0.707106781186547 ; |
309 |
|
|
fudge = 1./sqrt(2. * ((double) n ) ); |
310 |
|
|
for ( j = 0 ; j < n ; j ++ ) |
311 |
|
|
cos_data[j] *= fudge ; |
312 |
|
|
|
313 |
|
|
/* |
314 |
|
|
do the projections; Note that the cos_pml_table has |
315 |
|
|
had all the zeroes stripped out so the indexing is |
316 |
|
|
complicated somewhat |
317 |
|
|
*/ |
318 |
|
|
|
319 |
|
|
|
320 |
|
|
/******** this is the original loop |
321 |
|
|
|
322 |
|
|
toggle = 0 ; |
323 |
|
|
for (i=m; i<bw; i++) |
324 |
|
|
{ |
325 |
|
|
pml_ptr = cos_pml_table + NewTableOffset(m,i); |
326 |
|
|
|
327 |
|
|
if ((m % 2) == 0) |
328 |
|
|
{ |
329 |
|
|
for (j=0; j<(i/2)+1; j++) |
330 |
|
|
result[i-m] += cos_data[(2*j)+toggle] * pml_ptr[j]; |
331 |
|
|
} |
332 |
|
|
else |
333 |
|
|
{ |
334 |
|
|
if (((i-m) % 2) == 0) |
335 |
|
|
{ |
336 |
|
|
for (j=0; j<(i/2)+1; j++) |
337 |
|
|
result[i-m] += cos_data[(2*j)+toggle] * pml_ptr[j]; |
338 |
|
|
} |
339 |
|
|
else |
340 |
|
|
{ |
341 |
|
|
for (j=0; j<(i/2); j++) |
342 |
|
|
result[i-m] += cos_data[(2*j)+toggle] * pml_ptr[j]; |
343 |
|
|
} |
344 |
|
|
} |
345 |
|
|
|
346 |
|
|
toggle = (toggle+1) % 2; |
347 |
|
|
} |
348 |
|
|
|
349 |
|
|
*****/ |
350 |
|
|
|
351 |
|
|
/******** this is the new loop *********/ |
352 |
|
|
toggle = 0 ; |
353 |
|
|
for ( i=m; i<bw; i++ ) |
354 |
|
|
{ |
355 |
|
|
pml_ptr = cos_pml_table + NewTableOffset(m,i); |
356 |
|
|
|
357 |
|
|
result0 = 0.0 ; result1 = 0.0 ; |
358 |
|
|
result2 = 0.0 ; result3 = 0.0 ; |
359 |
|
|
|
360 |
|
|
for ( j = 0 ; j < ( (i/2) % 4 ) ; ++j ) |
361 |
|
|
result0 += cos_data[(2*j)+toggle] * pml_ptr[j]; |
362 |
|
|
|
363 |
|
|
for ( ; j < (i/2) ; j += 4 ) |
364 |
|
|
{ |
365 |
|
|
result0 += cos_data[(2*j)+toggle] * pml_ptr[j]; |
366 |
|
|
result1 += cos_data[(2*(j+1))+toggle] * pml_ptr[j+1]; |
367 |
|
|
result2 += cos_data[(2*(j+2))+toggle] * pml_ptr[j+2]; |
368 |
|
|
result3 += cos_data[(2*(j+3))+toggle] * pml_ptr[j+3]; |
369 |
|
|
} |
370 |
|
|
|
371 |
|
|
if ((((i-m) % 2) == 0 ) || ( (m % 2) == 0 )) |
372 |
|
|
result0 += cos_data[(2*(i/2))+toggle] * pml_ptr[(i/2)]; |
373 |
|
|
|
374 |
|
|
result[i-m] = result0 + result1 + result2 + result3 ; |
375 |
|
|
|
376 |
|
|
toggle = (toggle + 1)%2 ; |
377 |
|
|
|
378 |
|
|
} |
379 |
|
|
} |
380 |
|
|
|
381 |
|
|
|