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