00001 #include <stdio.h>
00002 #include <math.h>
00003 #include <string.h>
00004 #include <stdlib.h>
00005
00006 #include "util.h"
00007 #include "nfft3.h"
00008 #include "options.h"
00009
00010
00011
00012
00013 void short_nfft_trafo_2d(nfft_plan* ths, nfft_plan* plan_1d)
00014 {
00015 int j,k0;
00016 double omega;
00017
00018 for(j=0;j<ths->M_total;j++)
00019 {
00020 ths->f[j]= 0;
00021 plan_1d->x[j] = ths->x[ths->d * j + 1];
00022 }
00023
00024 for(k0=0;k0<ths->N[0];k0++)
00025 {
00026 plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
00027
00028 nfft_trafo(plan_1d);
00029
00030 for(j=0;j<ths->M_total;j++)
00031 {
00032 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
00033 ths->f[j] += plan_1d->f[j] * cexp( - I*2*PI*omega);
00034 }
00035 }
00036 }
00037
00038 void short_nfft_adjoint_2d(nfft_plan* ths, nfft_plan* plan_1d)
00039 {
00040 int j,k0;
00041 double omega;
00042
00043 for(j=0;j<ths->M_total;j++)
00044 plan_1d->x[j] = ths->x[ths->d * j + 1];
00045
00046 for(k0=0;k0<ths->N[0];k0++)
00047 {
00048 for(j=0;j<ths->M_total;j++)
00049 {
00050 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
00051 plan_1d->f[j] = ths->f[j] * cexp( + I*2*PI*omega);
00052 }
00053
00054 plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
00055
00056 nfft_adjoint(plan_1d);
00057 }
00058 }
00059
00060
00061
00062
00063 void short_nfft_trafo_3d_1(nfft_plan* ths, nfft_plan* plan_1d)
00064 {
00065 int j,k0,k1;
00066 double omega;
00067
00068 for(j=0;j<ths->M_total;j++)
00069 {
00070 ths->f[j] = 0;
00071 plan_1d->x[j] = ths->x[ths->d * j + 2];
00072 }
00073
00074 for(k0=0;k0<ths->N[0];k0++)
00075 for(k1=0;k1<ths->N[1];k1++)
00076 {
00077 plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
00078
00079 nfft_trafo(plan_1d);
00080
00081 for(j=0;j<ths->M_total;j++)
00082 {
00083 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0]
00084 + ((double)(k1 - ths->N[1]/2)) * ths->x[ths->d * j + 1];
00085 ths->f[j] += plan_1d->f[j] * cexp( - I*2*PI*omega);
00086 }
00087 }
00088 }
00089
00090 void short_nfft_adjoint_3d_1(nfft_plan* ths, nfft_plan* plan_1d)
00091 {
00092 int j,k0,k1;
00093 double omega;
00094
00095 for(j=0;j<ths->M_total;j++)
00096 plan_1d->x[j] = ths->x[ths->d * j + 2];
00097
00098 for(k0=0;k0<ths->N[0];k0++)
00099 for(k1=0;k1<ths->N[1];k1++)
00100 {
00101 for(j=0;j<ths->M_total;j++)
00102 {
00103 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0]
00104 + ((double)(k1 - ths->N[1]/2)) * ths->x[ths->d * j + 1];
00105 plan_1d->f[j] = ths->f[j] * cexp( + I*2*PI*omega);
00106 }
00107
00108 plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
00109
00110 nfft_adjoint(plan_1d);
00111 }
00112 }
00113
00114
00115
00116
00117 void short_nfft_trafo_3d_2(nfft_plan* ths, nfft_plan* plan_2d)
00118 {
00119 int j,k0;
00120 double omega;
00121
00122 for(j=0;j<ths->M_total;j++)
00123 {
00124 ths->f[j] = 0;
00125 plan_2d->x[2*j+0] = ths->x[ths->d * j + 1];
00126 plan_2d->x[2*j+1] = ths->x[ths->d * j + 2];
00127 }
00128
00129 for(k0=0;k0<ths->N[0];k0++)
00130 {
00131 plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
00132
00133 nfft_trafo(plan_2d);
00134
00135 for(j=0;j<ths->M_total;j++)
00136 {
00137 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
00138 ths->f[j] += plan_2d->f[j] * cexp( - I*2*PI*omega);
00139 }
00140 }
00141 }
00142
00143 void short_nfft_adjoint_3d_2(nfft_plan* ths, nfft_plan* plan_2d)
00144 {
00145 int j,k0;
00146 double omega;
00147
00148 for(j=0;j<ths->M_total;j++)
00149 {
00150 plan_2d->x[2*j+0] = ths->x[ths->d * j + 1];
00151 plan_2d->x[2*j+1] = ths->x[ths->d * j + 2];
00152 }
00153
00154 for(k0=0;k0<ths->N[0];k0++)
00155 {
00156 for(j=0;j<ths->M_total;j++)
00157 {
00158 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
00159 plan_2d->f[j] = ths->f[j] * cexp( + I*2*PI*omega);
00160 }
00161
00162 plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
00163
00164 nfft_adjoint(plan_2d);
00165 }
00166 }
00167
00168
00169
00170 int index_sparse_to_full_direct_2d(int J, int k)
00171 {
00172 int N=nfft_int_2_pow(J+2);
00173 int N_B=nfft_int_2_pow(J);
00174
00175 int j=k/N_B;
00176 int r=j/4;
00177
00178 int i, o, a, b,s,l,m1,m2;
00179 int k1,k2;
00180
00181 if (k>=(J+4)*nfft_int_2_pow(J+1))
00182 {
00183 printf("Fehler!\n");
00184 return(-1);
00185 }
00186 else
00187 {
00188 if (r>(J+1)/2)
00189 {
00190 i=k-4*((J+1)/2+1)*N_B;
00191 a=nfft_int_2_pow(J/2+1);
00192 m1=i/a;
00193 m2=i%a;
00194 k1=N/2-a/2+m1;
00195 k2=N/2-a/2+m2;
00196 }
00197 else
00198 {
00199 i=k-j*N_B;
00200 o=j%4;
00201 a=nfft_int_2_pow(r);
00202 b=nfft_int_2_pow(J-r);
00203 l=NFFT_MAX(a,b);
00204 s=NFFT_MIN(a,b);
00205 m1=i/l;
00206 m2=i%l;
00207
00208 switch(o)
00209 {
00210 case 0:
00211 {
00212 k1=N/2-a/2 ;
00213 k2=N/2+ b ;
00214
00215 if (b>=a)
00216 {
00217 k1+=m1;
00218 k2+=m2;
00219 }
00220 else
00221 {
00222 k1+=m2;
00223 k2+=m1;
00224 }
00225 break;
00226 }
00227 case 1:
00228 {
00229 k1=N/2+ b ;
00230 k2=N/2-a/2 ;
00231
00232 if (b>a)
00233 {
00234 k1+=m2;
00235 k2+=m1;
00236 }
00237 else
00238 {
00239 k1+=m1;
00240 k2+=m2;
00241 }
00242 break;
00243 }
00244 case 2:
00245 {
00246 k1=N/2-a/2 ;
00247 k2=N/2-2*b ;
00248
00249 if (b>=a)
00250 {
00251 k1+=m1;
00252 k2+=m2;
00253 }
00254 else
00255 {
00256 k1+=m2;
00257 k2+=m1;
00258 }
00259 break;
00260 }
00261 case 3:
00262 {
00263 k1=N/2-2*b ;
00264 k2=N/2-a/2 ;
00265
00266 if (b>a)
00267 {
00268 k1+=m2;
00269 k2+=m1;
00270 }
00271 else
00272 {
00273 k1+=m1;
00274 k2+=m2;
00275 }
00276 break;
00277 }
00278 default:
00279 {
00280 k1=-1;
00281 k2=-1;
00282 }
00283 }
00284 }
00285
00286 return(k1*N+k2);
00287 }
00288 }
00289
00290 inline int index_sparse_to_full_2d(nsfft_plan *ths, int k)
00291 {
00292
00293 if( k < ths->N_total)
00294 return ths->index_sparse_to_full[k];
00295 else
00296 return -1;
00297 }
00298
00299 int index_full_to_sparse_2d(int J, int k)
00300 {
00301 int N=nfft_int_2_pow(J+2);
00302 int N_B=nfft_int_2_pow(J);
00303
00304 int k1=k/N-N/2;
00305 int k2=k%N-N/2;
00306
00307 int r,a,b;
00308
00309 a=nfft_int_2_pow(J/2+1);
00310
00311 if ( (k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) )
00312 {
00313 return(4*((J+1)/2+1)*N_B+(k1+a/2)*a+(k2+a/2));
00314 }
00315
00316 for (r=0; r<=(J+1)/2; r++)
00317 {
00318 b=nfft_int_2_pow(r);
00319 a=nfft_int_2_pow(J-r);
00320 if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=a) && (k2<2*a) )
00321 {
00322 if (a>=b)
00323 return((4*r+0)*N_B+(k1+b/2)*a+(k2-a));
00324 else
00325 return((4*r+0)*N_B+(k2-a)*b+(k1+b/2));
00326 }
00327 else if ( (k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
00328 {
00329 if (a>b)
00330 return((4*r+1)*N_B+(k2+b/2)*a+(k1-a));
00331 else
00332 return((4*r+1)*N_B+(k1-a)*b+(k2+b/2));
00333 }
00334 else if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=-2*a) && (k2<-a) )
00335 {
00336 if (a>=b)
00337 return((4*r+2)*N_B+(k1+b/2)*a+(k2+2*a));
00338 else
00339 return((4*r+2)*N_B+(k2+2*a)*b+(k1+b/2));
00340 }
00341 else if ( (k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
00342 {
00343 if (a>b)
00344 return((4*r+3)*N_B+(k2+b/2)*a+(k1+2*a));
00345 else
00346 return((4*r+3)*N_B+(k1+2*a)*b+(k2+b/2));
00347 }
00348 }
00349
00350 return(-1);
00351 }
00352
00353 void init_index_sparse_to_full_2d(nsfft_plan *ths)
00354 {
00355 int k_S;
00356
00357 for (k_S=0; k_S<ths->N_total; k_S++)
00358 ths->index_sparse_to_full[k_S]=index_sparse_to_full_direct_2d(ths->J, k_S);
00359 }
00360
00361 inline int index_sparse_to_full_3d(nsfft_plan *ths, int k)
00362 {
00363
00364 if( k < ths->N_total)
00365 return ths->index_sparse_to_full[k];
00366 else
00367 return -1;
00368 }
00369
00370 int index_full_to_sparse_3d(int J, int k)
00371 {
00372 int N=nfft_int_2_pow(J+2);
00373 int N_B_r;
00374 int sum_N_B_less_r;
00375
00376 int r,a,b;
00377
00378 int k3=(k%N)-N/2;
00379 int k2=((k/N)%N)-N/2;
00380 int k1=k/(N*N)-N/2;
00381
00382 a=nfft_int_2_pow(J/2+1);
00383
00384 if((k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) && (k3>=(-a/2)) &&
00385 (k3<a/2))
00386 {
00387 return(6*nfft_int_2_pow(J)*(nfft_int_2_pow((J+1)/2+1)-1)+((k1+a/2)*a+(k2+a/2))*a+
00388 (k3+a/2));
00389 }
00390
00391 sum_N_B_less_r=0;
00392 for (r=0; r<=(J+1)/2; r++)
00393 {
00394 a=nfft_int_2_pow(J-r);
00395 b=nfft_int_2_pow(r);
00396
00397 N_B_r=a*b*b;
00398
00399
00400 if ((k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
00401 (k3>=-(b/2)) && (k3<(b+1)/2))
00402 {
00403 if(a>b)
00404 return sum_N_B_less_r+N_B_r*0 + ((k2+b/2)*b+k3+b/2)*a + (k1-a);
00405 else
00406 return sum_N_B_less_r+N_B_r*0 + ((k1-a)*b+(k2+b/2))*b + (k3+b/2);
00407 }
00408 else if ((k2>=a) && (k2<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
00409 (k3>=-(b/2)) && (k3<(b+1)/2))
00410 {
00411 if(a>b)
00412 return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+k3+b/2)*a + (k2-a);
00413 else if (a==b)
00414 return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+(k2-a))*a + (k3+b/2);
00415 else
00416 return sum_N_B_less_r+N_B_r*1 + ((k2-a)*b+(k1+b/2))*b + (k3+b/2);
00417 }
00418 else if ((k3>=a) && (k3<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
00419 (k2>=-(b/2)) && (k2<(b+1)/2))
00420 {
00421 if(a>=b)
00422 return sum_N_B_less_r+N_B_r*2 + ((k1+b/2)*b+k2+b/2)*a + (k3-a);
00423 else
00424 return sum_N_B_less_r+N_B_r*2 + ((k3-a)*b+(k1+b/2))*b + (k2+b/2);
00425 }
00426
00427 else if ((k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
00428 (k3>=-(b/2)) && (k3<(b+1)/2))
00429 {
00430 if(a>b)
00431 return sum_N_B_less_r+N_B_r*3 + ((k2+b/2)*b+k3+b/2)*a + (k1+2*a);
00432 else
00433 return sum_N_B_less_r+N_B_r*3 + ((k1+2*a)*b+(k2+b/2))*b + (k3+b/2);
00434 }
00435 else if ((k2>=-2*a) && (k2<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
00436 (k3>=-(b/2)) && (k3<(b+1)/2))
00437 {
00438 if(a>b)
00439 return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+k3+b/2)*a + (k2+2*a);
00440 else if (a==b)
00441 return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+(k2+2*a))*a + (k3+b/2);
00442 else
00443 return sum_N_B_less_r+N_B_r*4 + ((k2+2*a)*b+(k1+b/2))*b + (k3+b/2);
00444 }
00445 else if ((k3>=-2*a) && (k3<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
00446 (k2>=-(b/2)) && (k2<(b+1)/2))
00447 {
00448 if(a>=b)
00449 return sum_N_B_less_r+N_B_r*5 + ((k1+b/2)*b+k2+b/2)*a + (k3+2*a);
00450 else
00451 return sum_N_B_less_r+N_B_r*5 + ((k3+2*a)*b+(k1+b/2))*b + (k2+b/2);
00452 }
00453
00454 sum_N_B_less_r+=6*N_B_r;
00455 }
00456
00457 return(-1);
00458 }
00459
00460 void init_index_sparse_to_full_3d(nsfft_plan *ths)
00461 {
00462 int k1,k2,k3,k_s,r;
00463 int a,b;
00464 int N=nfft_int_2_pow(ths->J+2);
00465 int Nc=ths->center_nfft_plan->N[0];
00466
00467 for (k_s=0, r=0; r<=(ths->J+1)/2; r++)
00468 {
00469 a=nfft_int_2_pow(ths->J-r);
00470 b=nfft_int_2_pow(r);
00471
00472
00473
00474
00475 if(a>b)
00476 for(k2=-b/2;k2<(b+1)/2;k2++)
00477 for(k3=-b/2;k3<(b+1)/2;k3++)
00478 for(k1=a; k1<2*a; k1++,k_s++)
00479 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00480 else
00481 for(k1=a; k1<2*a; k1++)
00482 for(k2=-b/2;k2<(b+1)/2;k2++)
00483 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00484 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00485
00486
00487 if(a>b)
00488 for(k1=-b/2;k1<(b+1)/2;k1++)
00489 for(k3=-b/2;k3<(b+1)/2;k3++)
00490 for(k2=a; k2<2*a; k2++,k_s++)
00491 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00492 else if(a==b)
00493 for(k1=-b/2;k1<(b+1)/2;k1++)
00494 for(k2=a; k2<2*a; k2++)
00495 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00496 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00497 else
00498 for(k2=a; k2<2*a; k2++)
00499 for(k1=-b/2;k1<(b+1)/2;k1++)
00500 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00501 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00502
00503
00504 if(a>=b)
00505 for(k1=-b/2;k1<(b+1)/2;k1++)
00506 for(k2=-b/2;k2<(b+1)/2;k2++)
00507 for(k3=a; k3<2*a; k3++,k_s++)
00508 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00509 else
00510 for(k3=a; k3<2*a; k3++)
00511 for(k1=-b/2;k1<(b+1)/2;k1++)
00512 for(k2=-b/2;k2<(b+1)/2;k2++,k_s++)
00513 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00514
00515
00516 if(a>b)
00517 for(k2=-b/2;k2<(b+1)/2;k2++)
00518 for(k3=-b/2;k3<(b+1)/2;k3++)
00519 for(k1=-2*a; k1<-a; k1++,k_s++)
00520 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00521 else
00522 for(k1=-2*a; k1<-a; k1++)
00523 for(k2=-b/2;k2<(b+1)/2;k2++)
00524 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00525 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00526
00527
00528 if(a>b)
00529 for(k1=-b/2;k1<(b+1)/2;k1++)
00530 for(k3=-b/2;k3<(b+1)/2;k3++)
00531 for(k2=-2*a; k2<-a; k2++,k_s++)
00532 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00533 else if(a==b)
00534 for(k1=-b/2;k1<(b+1)/2;k1++)
00535 for(k2=-2*a; k2<-a; k2++)
00536 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00537 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00538 else
00539 for(k2=-2*a; k2<-a; k2++)
00540 for(k1=-b/2;k1<(b+1)/2;k1++)
00541 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
00542 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00543
00544
00545 if(a>=b)
00546 for(k1=-b/2;k1<(b+1)/2;k1++)
00547 for(k2=-b/2;k2<(b+1)/2;k2++)
00548 for(k3=-2*a; k3<-a; k3++,k_s++)
00549 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00550 else
00551 for(k3=-2*a; k3<-a; k3++)
00552 for(k1=-b/2;k1<(b+1)/2;k1++)
00553 for(k2=-b/2;k2<(b+1)/2;k2++,k_s++)
00554 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00555 }
00556
00557
00558 for(k1=-Nc/2;k1<Nc/2;k1++)
00559 for(k2=-Nc/2;k2<Nc/2;k2++)
00560 for(k3=-Nc/2; k3<Nc/2; k3++,k_s++)
00561 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
00562 }
00563
00564
00565 void nsfft_cp(nsfft_plan *ths, nfft_plan *ths_full_plan)
00566 {
00567 int k;
00568
00569
00570 memset(ths_full_plan->f_hat, 0, ths_full_plan->N_total*sizeof(double complex));
00571
00572
00573 for(k=0;k<ths->N_total;k++)
00574 ths_full_plan->f_hat[ths->index_sparse_to_full[k]]=ths->f_hat[k];
00575
00576
00577 memcpy(ths_full_plan->x,ths->act_nfft_plan->x,ths->M_total*ths->d*sizeof(double));
00578 }
00579
00580
00581 void test_copy_sparse_to_full_2d(nsfft_plan *ths, nfft_plan *ths_full_plan)
00582 {
00583 int r;
00584 int k1, k2;
00585 int a,b;
00586 const int J=ths->J;
00587 const int N=ths_full_plan->N[0];
00588 const int N_B=nfft_int_2_pow(J);
00589
00590
00591 nsfft_cp(ths, ths_full_plan);
00592
00593
00594 printf("f_hat blockwise\n");
00595 for (r=0; r<=(J+1)/2; r++){
00596 a=nfft_int_2_pow(J-r); b=nfft_int_2_pow(r);
00597
00598 printf("top\n");
00599 for (k1=0; k1<a; k1++){
00600 for (k2=0; k2<b; k2++){
00601 printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+1)*N_B+ k1*b+k2]),
00602 cimag(ths->f_hat[(4*r+1)*N_B+ k1*b+k2]));
00603 }
00604 printf("\n");
00605 }
00606
00607 printf("bottom\n");
00608 for (k1=0; k1<a; k1++){
00609 for (k2=0; k2<b; k2++){
00610 printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+3)*N_B+ k1*b+k2]),
00611 cimag(ths->f_hat[(4*r+3)*N_B+ k1*b+k2]));
00612 }
00613 printf("\n");
00614 }
00615
00616 printf("right\n");
00617 for (k2=0; k2<b; k2++){
00618 for (k1=0; k1<a; k1++){
00619 printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+0)*N_B+ k2*a+k1]),
00620 cimag(ths->f_hat[(4*r+0)*N_B+ k2*a+k1]));
00621 }
00622 printf("\n");
00623 }
00624
00625 printf("left\n");
00626 for (k2=0; k2<b; k2++){
00627 for (k1=0; k1<a; k1++){
00628 printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+2)*N_B+ k2*a+k1]),
00629 cimag(ths->f_hat[(4*r+2)*N_B+ k2*a+k1]));
00630 }
00631 printf("\n");
00632 }
00633 }
00634
00635 return;
00636
00637 printf("full f_hat\n");
00638 for (k1=0;k1<N;k1++){
00639 for (k2=0;k2<N;k2++){
00640 printf("(%1.1f,%1.1f) ", creal(ths_full_plan->f_hat[k1*N+k2]),
00641 cimag(ths_full_plan->f_hat[k1*N+k2]));
00642 }
00643 printf("\n");
00644 }
00645 }
00646
00647 void test_sparse_to_full_2d(nsfft_plan* ths)
00648 {
00649 int k_S,k1,k2;
00650 int N=nfft_int_2_pow(ths->J+2);
00651
00652 printf("N=%d\n\n",N);
00653
00654 for(k1=0;k1<N;k1++)
00655 for(k2=0;k2<N;k2++)
00656 {
00657 k_S=index_full_to_sparse_2d(ths->J, k1*N+k2);
00658 if(k_S!=-1)
00659 printf("(%+d, %+d)\t= %+d \t= %+d = %+d \n",k1-N/2,k2-N/2,
00660 k1*N+k2, k_S, ths->index_sparse_to_full[k_S]);
00661 }
00662 }
00663
00664 void test_sparse_to_full_3d(nsfft_plan* ths)
00665 {
00666 int k_S,k1,k2,k3;
00667 int N=nfft_int_2_pow(ths->J+2);
00668
00669 printf("N=%d\n\n",N);
00670
00671 for(k1=0;k1<N;k1++)
00672 for(k2=0;k2<N;k2++)
00673 for(k3=0;k3<N;k3++)
00674 {
00675 k_S=index_full_to_sparse_3d(ths->J, (k1*N+k2)*N+k3);
00676 if(k_S!=-1)
00677 printf("(%d, %d, %d)\t= %d \t= %d = %d \n",k1-N/2,k2-N/2,k3-N/2,
00678 (k1*N+k2)*N+k3,k_S, ths->index_sparse_to_full[k_S]);
00679 }
00680 }
00681
00682
00683 void nsfft_init_random_nodes_coeffs(nsfft_plan *ths)
00684 {
00685 int j;
00686
00687
00688 nfft_vrand_unit_complex(ths->f_hat, ths->N_total);
00689
00690
00691 nfft_vrand_shifted_unit_double(ths->act_nfft_plan->x, ths->d * ths->M_total);
00692
00693 if(ths->d==2)
00694 for(j=0;j<ths->M_total;j++)
00695 {
00696 ths->x_transposed[2*j+0]=ths->act_nfft_plan->x[2*j+1];
00697 ths->x_transposed[2*j+1]=ths->act_nfft_plan->x[2*j+0];
00698 }
00699 else
00700 for(j=0;j<ths->M_total;j++)
00701 {
00702 ths->x_102[3*j+0]=ths->act_nfft_plan->x[3*j+1];
00703 ths->x_102[3*j+1]=ths->act_nfft_plan->x[3*j+0];
00704 ths->x_102[3*j+2]=ths->act_nfft_plan->x[3*j+2];
00705
00706 ths->x_201[3*j+0]=ths->act_nfft_plan->x[3*j+2];
00707 ths->x_201[3*j+1]=ths->act_nfft_plan->x[3*j+0];
00708 ths->x_201[3*j+2]=ths->act_nfft_plan->x[3*j+1];
00709
00710 ths->x_120[3*j+0]=ths->act_nfft_plan->x[3*j+1];
00711 ths->x_120[3*j+1]=ths->act_nfft_plan->x[3*j+2];
00712 ths->x_120[3*j+2]=ths->act_nfft_plan->x[3*j+0];
00713
00714 ths->x_021[3*j+0]=ths->act_nfft_plan->x[3*j+0];
00715 ths->x_021[3*j+1]=ths->act_nfft_plan->x[3*j+2];
00716 ths->x_021[3*j+2]=ths->act_nfft_plan->x[3*j+1];
00717 }
00718 }
00719
00720 void nsdft_trafo_2d(nsfft_plan *ths)
00721 {
00722 int j,k_S,k_L,k0,k1;
00723 double omega;
00724 int N=nfft_int_2_pow(ths->J+2);
00725
00726 memset(ths->f,0,ths->M_total*sizeof(double complex));
00727
00728 for(k_S=0;k_S<ths->N_total;k_S++)
00729 {
00730 k_L=ths->index_sparse_to_full[k_S];
00731 k0=k_L / N;
00732 k1=k_L % N;
00733
00734 for(j=0;j<ths->M_total;j++)
00735 {
00736 omega =
00737 ((double)(k0 - N/2)) * ths->act_nfft_plan->x[2 * j + 0] +
00738 ((double)(k1 - N/2)) * ths->act_nfft_plan->x[2 * j + 1];
00739 ths->f[j] += ths->f_hat[k_S] * cexp( - I*2*PI*omega);
00740 }
00741 }
00742 }
00743
00744 void nsdft_trafo_3d(nsfft_plan *ths)
00745 {
00746 int j,k_S,k0,k1,k2;
00747 double omega;
00748 int N=nfft_int_2_pow(ths->J+2);
00749 int k_L;
00750
00751 memset(ths->f,0,ths->M_total*sizeof(double complex));
00752
00753 for(k_S=0;k_S<ths->N_total;k_S++)
00754 {
00755 k_L=ths->index_sparse_to_full[k_S];
00756
00757 k0=k_L/(N*N);
00758 k1=(k_L/N)%N;
00759 k2=k_L%N;
00760
00761 for(j=0;j<ths->M_total;j++)
00762 {
00763 omega =
00764 ((double)(k0 - N/2)) * ths->act_nfft_plan->x[3 * j + 0] +
00765 ((double)(k1 - N/2)) * ths->act_nfft_plan->x[3 * j + 1] +
00766 ((double)(k2 - N/2)) * ths->act_nfft_plan->x[3 * j + 2];
00767 ths->f[j] += ths->f_hat[k_S] * cexp( - I*2*PI*omega);
00768 }
00769 }
00770 }
00771
00772 void nsdft_trafo(nsfft_plan *ths)
00773 {
00774 if(ths->d==2)
00775 nsdft_trafo_2d(ths);
00776 else
00777 nsdft_trafo_3d(ths);
00778 }
00779
00780 void nsdft_adjoint_2d(nsfft_plan *ths)
00781 {
00782 int j,k_S,k_L,k0,k1;
00783 double omega;
00784 int N=nfft_int_2_pow(ths->J+2);
00785
00786 memset(ths->f_hat,0,ths->N_total*sizeof(double complex));
00787
00788 for(k_S=0;k_S<ths->N_total;k_S++)
00789 {
00790 k_L=ths->index_sparse_to_full[k_S];
00791 k0=k_L / N;
00792 k1=k_L % N;
00793
00794 for(j=0;j<ths->M_total;j++)
00795 {
00796 omega =
00797 ((double)(k0 - N/2)) * ths->act_nfft_plan->x[2 * j + 0] +
00798 ((double)(k1 - N/2)) * ths->act_nfft_plan->x[2 * j + 1];
00799 ths->f_hat[k_S] += ths->f[j] * cexp( + I*2*PI*omega);
00800 }
00801 }
00802 }
00803
00804 void nsdft_adjoint_3d(nsfft_plan *ths)
00805 {
00806 int j,k_S,k0,k1,k2;
00807 double omega;
00808 int N=nfft_int_2_pow(ths->J+2);
00809 int k_L;
00810
00811 memset(ths->f_hat,0,ths->N_total*sizeof(double complex));
00812
00813 for(k_S=0;k_S<ths->N_total;k_S++)
00814 {
00815 k_L=ths->index_sparse_to_full[k_S];
00816
00817 k0=k_L/(N*N);
00818 k1=(k_L/N)%N;
00819 k2=k_L%N;
00820
00821 for(j=0;j<ths->M_total;j++)
00822 {
00823 omega =
00824 ((double)(k0 - N/2)) * ths->act_nfft_plan->x[3 * j + 0] +
00825 ((double)(k1 - N/2)) * ths->act_nfft_plan->x[3 * j + 1] +
00826 ((double)(k2 - N/2)) * ths->act_nfft_plan->x[3 * j + 2];
00827 ths->f_hat[k_S] += ths->f[j] * cexp( + I*2*PI*omega);
00828 }
00829 }
00830 }
00831
00832 void nsdft_adjoint(nsfft_plan *ths)
00833 {
00834 if(ths->d==2)
00835 nsdft_adjoint_2d(ths);
00836 else
00837 nsdft_adjoint_3d(ths);
00838 }
00839
00840 void nsfft_trafo_2d(nsfft_plan *ths)
00841 {
00842 int r,rr,j;
00843 double temp;
00844
00845 int M=ths->M_total;
00846 int J=ths->J;
00847
00848
00849 ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*nfft_int_2_pow(J);
00850
00851 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
00852 ndft_trafo(ths->center_nfft_plan);
00853 else
00854 nfft_trafo(ths->center_nfft_plan);
00855
00856 for (j=0; j<M; j++)
00857 ths->f[j] = ths->center_nfft_plan->f[j];
00858
00859 for(rr=0;rr<=(J+1)/2;rr++)
00860 {
00861 r=NFFT_MIN(rr,J-rr);
00862 ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[r];
00863 ths->act_nfft_plan->N[0]=nfft_int_2_pow(r); ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
00864 ths->act_nfft_plan->N[1]=nfft_int_2_pow(J-r); ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
00865
00866
00867
00868 temp=-3.0*PI*nfft_int_2_pow(J-rr);
00869
00870
00871 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*nfft_int_2_pow(J);
00872
00873 if(r<rr)
00874 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00875
00876 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00877 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00878 ndft_trafo(ths->act_nfft_plan);
00879 else
00880 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
00881 else
00882 nfft_trafo(ths->act_nfft_plan);
00883
00884 if(r<rr)
00885 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00886
00887 for (j=0; j<M; j++)
00888 ths->f[j] += ths->act_nfft_plan->f[j] *
00889 cexp( + I*temp*ths->act_nfft_plan->x[2*j+1]);
00890
00891
00892 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*nfft_int_2_pow(J);
00893
00894 if((r==rr)&&(J-rr!=rr))
00895 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00896
00897 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00898 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00899 ndft_trafo(ths->act_nfft_plan);
00900 else
00901 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
00902 else
00903 nfft_trafo(ths->act_nfft_plan);
00904
00905 if((r==rr)&&(J-rr!=rr))
00906 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00907
00908 for (j=0; j<M; j++)
00909 ths->f[j] += ths->act_nfft_plan->f[j] *
00910 cexp( + I*temp*ths->act_nfft_plan->x[2*j+0]);
00911
00912
00913 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*nfft_int_2_pow(J);
00914
00915 if(r<rr)
00916 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00917
00918 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00919 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00920 ndft_trafo(ths->act_nfft_plan);
00921 else
00922 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
00923 else
00924 nfft_trafo(ths->act_nfft_plan);
00925
00926 if(r<rr)
00927 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00928
00929 for (j=0; j<M; j++)
00930 ths->f[j] += ths->act_nfft_plan->f[j] *
00931 cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
00932
00933
00934 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*nfft_int_2_pow(J);
00935
00936 if((r==rr)&&(J-rr!=rr))
00937 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00938
00939 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00940 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00941 ndft_trafo(ths->act_nfft_plan);
00942 else
00943 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
00944 else
00945 nfft_trafo(ths->act_nfft_plan);
00946
00947 if((r==rr)&&(J-rr!=rr))
00948 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00949
00950 for (j=0; j<M; j++)
00951 ths->f[j] += ths->act_nfft_plan->f[j] *
00952 cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
00953 }
00954 }
00955
00956 void nsfft_adjoint_2d(nsfft_plan *ths)
00957 {
00958 int r,rr,j;
00959 double temp;
00960
00961 int M=ths->M_total;
00962 int J=ths->J;
00963
00964
00965 for (j=0; j<M; j++)
00966 ths->center_nfft_plan->f[j] = ths->f[j];
00967
00968 ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*nfft_int_2_pow(J);
00969
00970 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
00971 ndft_adjoint(ths->center_nfft_plan);
00972 else
00973 nfft_adjoint(ths->center_nfft_plan);
00974
00975 for(rr=0;rr<=(J+1)/2;rr++)
00976 {
00977 r=NFFT_MIN(rr,J-rr);
00978 ths->act_nfft_plan->my_fftw_plan2 = ths->set_fftw_plan2[r];
00979 ths->act_nfft_plan->N[0]=nfft_int_2_pow(r); ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
00980 ths->act_nfft_plan->N[1]=nfft_int_2_pow(J-r); ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
00981
00982
00983
00984 temp=-3.0*PI*nfft_int_2_pow(J-rr);
00985
00986
00987 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*nfft_int_2_pow(J);
00988
00989 for (j=0; j<M; j++)
00990 ths->act_nfft_plan->f[j]= ths->f[j] *
00991 cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
00992
00993 if(r<rr)
00994 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
00995
00996 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
00997 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
00998 ndft_adjoint(ths->act_nfft_plan);
00999 else
01000 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01001 else
01002 nfft_adjoint(ths->act_nfft_plan);
01003
01004 if(r<rr)
01005 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
01006
01007
01008 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*nfft_int_2_pow(J);
01009
01010 for (j=0; j<M; j++)
01011 ths->act_nfft_plan->f[j]= ths->f[j] *
01012 cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
01013
01014 if((r==rr)&&(J-rr!=rr))
01015 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
01016
01017 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01018 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01019 ndft_adjoint(ths->act_nfft_plan);
01020 else
01021 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01022 else
01023 nfft_adjoint(ths->act_nfft_plan);
01024
01025 if((r==rr)&&(J-rr!=rr))
01026 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
01027
01028
01029 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*nfft_int_2_pow(J);
01030
01031 for (j=0; j<M; j++)
01032 ths->act_nfft_plan->f[j]= ths->f[j] *
01033 cexp( + I*temp*ths->act_nfft_plan->x[2*j+1]);
01034
01035 if(r<rr)
01036 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
01037
01038 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01039 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01040 ndft_adjoint(ths->act_nfft_plan);
01041 else
01042 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01043 else
01044 nfft_adjoint(ths->act_nfft_plan);
01045
01046 if(r<rr)
01047 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
01048
01049
01050 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*nfft_int_2_pow(J);
01051
01052 for (j=0; j<M; j++)
01053 ths->act_nfft_plan->f[j]= ths->f[j] *
01054 cexp( + I*temp*ths->act_nfft_plan->x[2*j+0]);
01055
01056 if((r==rr)&&(J-rr!=rr))
01057 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
01058
01059 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01060 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01061 ndft_adjoint(ths->act_nfft_plan);
01062 else
01063 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01064 else
01065 nfft_adjoint(ths->act_nfft_plan);
01066
01067 if((r==rr)&&(J-rr!=rr))
01068 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_transposed);
01069 }
01070 }
01071
01072 void nsfft_trafo_3d(nsfft_plan *ths)
01073 {
01074 int r,rr,j;
01075 double temp;
01076 int sum_N_B_less_r,N_B_r,a,b;
01077
01078 int M=ths->M_total;
01079 int J=ths->J;
01080
01081
01082 ths->center_nfft_plan->f_hat=ths->f_hat+6*nfft_int_2_pow(J)*(nfft_int_2_pow((J+1)/2+1)-1);
01083
01084 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
01085 ndft_trafo(ths->center_nfft_plan);
01086 else
01087 nfft_trafo(ths->center_nfft_plan);
01088
01089 for (j=0; j<M; j++)
01090 ths->f[j] = ths->center_nfft_plan->f[j];
01091
01092 sum_N_B_less_r=0;
01093 for(rr=0;rr<=(J+1)/2;rr++)
01094 {
01095 a=nfft_int_2_pow(J-rr);
01096 b=nfft_int_2_pow(rr);
01097
01098 N_B_r=a*b*b;
01099
01100 r=NFFT_MIN(rr,J-rr);
01101 ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
01102
01103 ths->act_nfft_plan->N[0]=nfft_int_2_pow(r);
01104 if(a<b)
01105 ths->act_nfft_plan->N[1]=nfft_int_2_pow(J-r);
01106 else
01107 ths->act_nfft_plan->N[1]=nfft_int_2_pow(r);
01108 ths->act_nfft_plan->N[2]=nfft_int_2_pow(J-r);
01109
01110
01111
01112 ths->act_nfft_plan->N_total=ths->act_nfft_plan->N[0]*ths->act_nfft_plan->N[1]*ths->act_nfft_plan->N[2];
01113 ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
01114 ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
01115 ths->act_nfft_plan->n[2]=ths->sigma*ths->act_nfft_plan->N[2];
01116 ths->act_nfft_plan->n_total=ths->act_nfft_plan->n[0]*ths->act_nfft_plan->n[1]*ths->act_nfft_plan->n[2];
01117
01118
01119 if((J==0)||((J==1)&&(rr==1)))
01120 temp=-2.0*PI;
01121 else
01122 temp=-3.0*PI*nfft_int_2_pow(J-rr);
01123
01124
01125 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
01126
01127 if(a>b)
01128 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_120);
01129
01130 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01131 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01132 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01133 ndft_trafo(ths->act_nfft_plan);
01134 else
01135 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01136 else
01137 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01138 else
01139 nfft_trafo(ths->act_nfft_plan);
01140
01141 if(a>b)
01142 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_120);
01143
01144 for (j=0; j<M; j++)
01145 ths->f[j] += ths->act_nfft_plan->f[j] *
01146 cexp( + I*temp*ths->act_nfft_plan->x[3*j+0]);
01147
01148
01149 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
01150
01151 if(a>b)
01152 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_021);
01153 if(a<b)
01154 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_102);
01155
01156 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01157 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01158 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01159 ndft_trafo(ths->act_nfft_plan);
01160 else
01161 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01162 else
01163 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01164 else
01165 nfft_trafo(ths->act_nfft_plan);
01166
01167 if(a>b)
01168 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_021);
01169 if(a<b)
01170 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_102);
01171
01172 for (j=0; j<M; j++)
01173 ths->f[j] += ths->act_nfft_plan->f[j] *
01174 cexp( + I*temp*ths->act_nfft_plan->x[3*j+1]);
01175
01176
01177 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
01178
01179 if(a<b)
01180 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_201);
01181
01182 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01183 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01184 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01185 ndft_trafo(ths->act_nfft_plan);
01186 else
01187 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01188 else
01189 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01190 else
01191 nfft_trafo(ths->act_nfft_plan);
01192
01193 if(a<b)
01194 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_201);
01195
01196 for (j=0; j<M; j++)
01197 ths->f[j] += ths->act_nfft_plan->f[j] *
01198 cexp( + I*temp*ths->act_nfft_plan->x[3*j+2]);
01199
01200
01201 if((J==0)||((J==1)&&(rr==1)))
01202 temp=-4.0*PI;
01203 else
01204 temp=-3.0*PI*nfft_int_2_pow(J-rr);
01205
01206
01207 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
01208
01209 if(a>b)
01210 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_120);
01211
01212 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01213 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01214 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01215 ndft_trafo(ths->act_nfft_plan);
01216 else
01217 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01218 else
01219 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01220 else
01221 nfft_trafo(ths->act_nfft_plan);
01222
01223 if(a>b)
01224 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_120);
01225
01226 for (j=0; j<M; j++)
01227 ths->f[j] += ths->act_nfft_plan->f[j] *
01228 cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
01229
01230
01231 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
01232
01233 if(a>b)
01234 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_021);
01235 if(a<b)
01236 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_102);
01237
01238 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01239 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01240 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01241 ndft_trafo(ths->act_nfft_plan);
01242 else
01243 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01244 else
01245 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01246 else
01247 nfft_trafo(ths->act_nfft_plan);
01248
01249 if(a>b)
01250 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_021);
01251 if(a<b)
01252 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_102);
01253
01254 for (j=0; j<M; j++)
01255 ths->f[j] += ths->act_nfft_plan->f[j] *
01256 cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
01257
01258
01259 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
01260
01261 if(a<b)
01262 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_201);
01263
01264 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01265 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01266 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01267 ndft_trafo(ths->act_nfft_plan);
01268 else
01269 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01270 else
01271 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01272 else
01273 nfft_trafo(ths->act_nfft_plan);
01274
01275 if(a<b)
01276 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_201);
01277
01278 for (j=0; j<M; j++)
01279 ths->f[j] += ths->act_nfft_plan->f[j] *
01280 cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
01281
01282 sum_N_B_less_r+=6*N_B_r;
01283 }
01284 }
01285
01286 void nsfft_adjoint_3d(nsfft_plan *ths)
01287 {
01288 int r,rr,j;
01289 double temp;
01290 int sum_N_B_less_r,N_B_r,a,b;
01291
01292 int M=ths->M_total;
01293 int J=ths->J;
01294
01295
01296 for (j=0; j<M; j++)
01297 ths->center_nfft_plan->f[j] = ths->f[j];
01298
01299 ths->center_nfft_plan->f_hat=ths->f_hat+6*nfft_int_2_pow(J)*(nfft_int_2_pow((J+1)/2+1)-1);
01300
01301 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
01302 ndft_adjoint(ths->center_nfft_plan);
01303 else
01304 nfft_adjoint(ths->center_nfft_plan);
01305
01306 sum_N_B_less_r=0;
01307 for(rr=0;rr<=(J+1)/2;rr++)
01308 {
01309 a=nfft_int_2_pow(J-rr);
01310 b=nfft_int_2_pow(rr);
01311
01312 N_B_r=a*b*b;
01313
01314 r=NFFT_MIN(rr,J-rr);
01315 ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
01316 ths->act_nfft_plan->my_fftw_plan2 = ths->set_fftw_plan2[rr];
01317
01318 ths->act_nfft_plan->N[0]=nfft_int_2_pow(r);
01319 if(a<b)
01320 ths->act_nfft_plan->N[1]=nfft_int_2_pow(J-r);
01321 else
01322 ths->act_nfft_plan->N[1]=nfft_int_2_pow(r);
01323 ths->act_nfft_plan->N[2]=nfft_int_2_pow(J-r);
01324
01325
01326
01327 ths->act_nfft_plan->N_total=ths->act_nfft_plan->N[0]*ths->act_nfft_plan->N[1]*ths->act_nfft_plan->N[2];
01328 ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
01329 ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
01330 ths->act_nfft_plan->n[2]=ths->sigma*ths->act_nfft_plan->N[2];
01331 ths->act_nfft_plan->n_total=ths->act_nfft_plan->n[0]*ths->act_nfft_plan->n[1]*ths->act_nfft_plan->n[2];
01332
01333
01334 if((J==0)||((J==1)&&(rr==1)))
01335 temp=-2.0*PI;
01336 else
01337 temp=-3.0*PI*nfft_int_2_pow(J-rr);
01338
01339
01340 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
01341
01342 for (j=0; j<M; j++)
01343 ths->act_nfft_plan->f[j]= ths->f[j] *
01344 cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
01345
01346 if(a>b)
01347 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_120);
01348
01349 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01350 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01351 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01352 ndft_adjoint(ths->act_nfft_plan);
01353 else
01354 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01355 else
01356 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01357 else
01358 nfft_adjoint(ths->act_nfft_plan);
01359
01360 if(a>b)
01361 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_120);
01362
01363
01364 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
01365
01366 for (j=0; j<M; j++)
01367 ths->act_nfft_plan->f[j]= ths->f[j] *
01368 cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
01369
01370 if(a>b)
01371 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_021);
01372 if(a<b)
01373 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_102);
01374
01375 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01376 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01377 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01378 ndft_adjoint(ths->act_nfft_plan);
01379 else
01380 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01381 else
01382 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01383 else
01384 nfft_adjoint(ths->act_nfft_plan);
01385
01386 if(a>b)
01387 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_021);
01388 if(a<b)
01389 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_102);
01390
01391
01392 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
01393
01394 for (j=0; j<M; j++)
01395 ths->act_nfft_plan->f[j]= ths->f[j] *
01396 cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
01397
01398 if(a<b)
01399 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_201);
01400
01401 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01402 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01403 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01404 ndft_adjoint(ths->act_nfft_plan);
01405 else
01406 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01407 else
01408 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01409 else
01410 nfft_adjoint(ths->act_nfft_plan);
01411
01412 if(a<b)
01413 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_201);
01414
01415
01416 if((J==0)||((J==1)&&(rr==1)))
01417 temp=-4.0*PI;
01418 else
01419 temp=-3.0*PI*nfft_int_2_pow(J-rr);
01420
01421
01422 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
01423
01424 for (j=0; j<M; j++)
01425 ths->act_nfft_plan->f[j]= ths->f[j] *
01426 cexp( + I*temp*ths->act_nfft_plan->x[3*j+0]);
01427
01428 if(a>b)
01429 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_120);
01430
01431 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01432 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01433 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01434 ndft_adjoint(ths->act_nfft_plan);
01435 else
01436 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01437 else
01438 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01439 else
01440 nfft_adjoint(ths->act_nfft_plan);
01441
01442 if(a>b)
01443 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_120);
01444
01445
01446 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
01447
01448 for (j=0; j<M; j++)
01449 ths->act_nfft_plan->f[j]= ths->f[j] *
01450 cexp( + I*temp*ths->act_nfft_plan->x[3*j+1]);
01451
01452 if(a>b)
01453 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_021);
01454 if(a<b)
01455 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_102);
01456
01457 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01458 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01459 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01460 ndft_adjoint(ths->act_nfft_plan);
01461 else
01462 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01463 else
01464 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01465 else
01466 nfft_adjoint(ths->act_nfft_plan);
01467
01468 if(a>b)
01469 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_021);
01470 if(a<b)
01471 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_102);
01472
01473
01474 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
01475
01476 for (j=0; j<M; j++)
01477 ths->act_nfft_plan->f[j]= ths->f[j] *
01478 cexp( + I*temp*ths->act_nfft_plan->x[3*j+2]);
01479
01480 if(a<b)
01481 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_201);
01482
01483 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
01484 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
01485 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
01486 ndft_adjoint(ths->act_nfft_plan);
01487 else
01488 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
01489 else
01490 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
01491 else
01492 nfft_adjoint(ths->act_nfft_plan);
01493
01494 if(a<b)
01495 NFFT_SWAP_double(ths->act_nfft_plan->x,ths->x_201);
01496
01497 sum_N_B_less_r+=6*N_B_r;
01498 }
01499 }
01500
01501 void nsfft_trafo(nsfft_plan *ths)
01502 {
01503 if(ths->d==2)
01504 nsfft_trafo_2d(ths);
01505 else
01506 nsfft_trafo_3d(ths);
01507 }
01508
01509 void nsfft_adjoint(nsfft_plan *ths)
01510 {
01511 if(ths->d==2)
01512 nsfft_adjoint_2d(ths);
01513 else
01514 nsfft_adjoint_3d(ths);
01515 }
01516
01517
01518
01519
01520 void nsfft_init_2d(nsfft_plan *ths, int J, int M, int m, unsigned snfft_flags)
01521 {
01522 int r;
01523 int N[2];
01524 int n[2];
01525
01526 ths->flags=snfft_flags;
01527 ths->sigma=2;
01528 ths->J=J;
01529 ths->M_total=M;
01530 ths->N_total=(J+4)*nfft_int_2_pow(J+1);
01531
01532
01533 ths->f = (double complex *)fftw_malloc(M*sizeof(double complex));
01534 ths->f_hat = (double complex *)fftw_malloc(ths->N_total*sizeof(double complex));
01535 ths->x_transposed= (double*)fftw_malloc(2*M*sizeof(double));
01536
01537 ths->act_nfft_plan = (nfft_plan*)fftw_malloc(sizeof(nfft_plan));
01538 ths->center_nfft_plan = (nfft_plan*)fftw_malloc(sizeof(nfft_plan));
01539
01540 ths->set_fftw_plan1=(fftw_plan*) fftw_malloc((J/2+1)*sizeof(fftw_plan));
01541 ths->set_fftw_plan2=(fftw_plan*) fftw_malloc((J/2+1)*sizeof(fftw_plan));
01542
01543 ths->set_nfft_plan_1d = (nfft_plan*) fftw_malloc((nfft_ld(m)+1)*(sizeof(nfft_plan)));
01544
01545
01546
01547 N[0]=1; n[0]=ths->sigma*N[0];
01548 N[1]=nfft_int_2_pow(J); n[1]=ths->sigma*N[1];
01549
01550 nfft_init_guru(ths->act_nfft_plan,2,N,M,n,m, FG_PSI| MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
01551
01552 if(ths->act_nfft_plan->nfft_flags & PRE_ONE_PSI)
01553 nfft_precompute_one_psi(ths->act_nfft_plan);
01554
01555 ths->set_fftw_plan1[0]=ths->act_nfft_plan->my_fftw_plan1;
01556 ths->set_fftw_plan2[0]=ths->act_nfft_plan->my_fftw_plan2;
01557
01558 for(r=1;r<=J/2;r++)
01559 {
01560 N[0]=nfft_int_2_pow(r); n[0]=ths->sigma*N[0];
01561 N[1]=nfft_int_2_pow(J-r); n[1]=ths->sigma*N[1];
01562 ths->set_fftw_plan1[r] =
01563 fftw_plan_dft(2, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
01564 FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
01565
01566 ths->set_fftw_plan2[r] =
01567 fftw_plan_dft(2, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
01568 FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
01569 }
01570
01571
01572 for(r=0;r<=nfft_ld(m);r++)
01573 {
01574 N[0]=nfft_int_2_pow(J-r); n[0]=ths->sigma*N[0];
01575
01576 nfft_init_guru(&(ths->set_nfft_plan_1d[r]),1,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
01577 ths->set_nfft_plan_1d[r].nfft_flags = ths->set_nfft_plan_1d[r].nfft_flags | FG_PSI;
01578 ths->set_nfft_plan_1d[r].K=ths->act_nfft_plan->K;
01579 ths->set_nfft_plan_1d[r].psi=ths->act_nfft_plan->psi;
01580 }
01581
01582
01583
01584 N[0]=nfft_int_2_pow(J/2+1); n[0]=ths->sigma*N[0];
01585 N[1]=nfft_int_2_pow(J/2+1); n[1]=ths->sigma*N[1];
01586 nfft_init_guru(ths->center_nfft_plan,2,N,M,n, m, MALLOC_F| FFTW_INIT,
01587 FFTW_MEASURE);
01588 ths->center_nfft_plan->x= ths->act_nfft_plan->x;
01589 ths->center_nfft_plan->nfft_flags = ths->center_nfft_plan->nfft_flags|
01590 FG_PSI;
01591 ths->center_nfft_plan->K=ths->act_nfft_plan->K;
01592 ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
01593
01594 if(ths->flags & NSDFT)
01595 {
01596 ths->index_sparse_to_full=(int*)fftw_malloc(ths->N_total*sizeof(int));
01597 init_index_sparse_to_full_2d(ths);
01598 }
01599 }
01600
01601
01602
01603 void nsfft_init_3d(nsfft_plan *ths, int J, int M, int m, unsigned snfft_flags)
01604 {
01605 int r,rr,a,b;
01606 int N[3];
01607 int n[3];
01608
01609 ths->flags=snfft_flags;
01610 ths->sigma=2;
01611 ths->J=J;
01612 ths->M_total=M;
01613 ths->N_total=6*nfft_int_2_pow(J)*(nfft_int_2_pow((J+1)/2+1)-1)+nfft_int_2_pow(3*(J/2+1));
01614
01615
01616 ths->f = (double complex *)fftw_malloc(M*sizeof(double complex));
01617 ths->f_hat = (double complex *)fftw_malloc(ths->N_total*sizeof(double complex));
01618
01619 ths->x_102= (double*)fftw_malloc(3*M*sizeof(double));
01620 ths->x_201= (double*)fftw_malloc(3*M*sizeof(double));
01621 ths->x_120= (double*)fftw_malloc(3*M*sizeof(double));
01622 ths->x_021= (double*)fftw_malloc(3*M*sizeof(double));
01623
01624 ths->act_nfft_plan = (nfft_plan*)fftw_malloc(sizeof(nfft_plan));
01625 ths->center_nfft_plan = (nfft_plan*)fftw_malloc(sizeof(nfft_plan));
01626
01627 ths->set_fftw_plan1=(fftw_plan*) fftw_malloc(((J+1)/2+1)*sizeof(fftw_plan));
01628 ths->set_fftw_plan2=(fftw_plan*) fftw_malloc(((J+1)/2+1)*sizeof(fftw_plan));
01629
01630 ths->set_nfft_plan_1d = (nfft_plan*) fftw_malloc((nfft_ld(m)+1)*(sizeof(nfft_plan)));
01631 ths->set_nfft_plan_2d = (nfft_plan*) fftw_malloc((nfft_ld(m)+1)*(sizeof(nfft_plan)));
01632
01633
01634
01635 N[0]=1; n[0]=ths->sigma*N[0];
01636 N[1]=1; n[1]=ths->sigma*N[1];
01637 N[2]=nfft_int_2_pow(J); n[2]=ths->sigma*N[2];
01638
01639 nfft_init_guru(ths->act_nfft_plan,3,N,M,n,m, FG_PSI| MALLOC_X| MALLOC_F, FFTW_MEASURE);
01640
01641 if(ths->act_nfft_plan->nfft_flags & PRE_ONE_PSI)
01642 nfft_precompute_one_psi(ths->act_nfft_plan);
01643
01644
01645 ths->act_nfft_plan->g1 = fftw_malloc(ths->sigma*ths->sigma*ths->sigma*nfft_int_2_pow(J+(J+1)/2)*sizeof(double complex));
01646 ths->act_nfft_plan->g2 = fftw_malloc(ths->sigma*ths->sigma*ths->sigma*nfft_int_2_pow(J+(J+1)/2)*sizeof(double complex));
01647
01648 ths->act_nfft_plan->my_fftw_plan1 =
01649 fftw_plan_dft(3, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
01650 FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
01651 ths->act_nfft_plan->my_fftw_plan2 =
01652 fftw_plan_dft(3, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
01653 FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
01654
01655 ths->set_fftw_plan1[0]=ths->act_nfft_plan->my_fftw_plan1;
01656 ths->set_fftw_plan2[0]=ths->act_nfft_plan->my_fftw_plan2;
01657
01658 for(rr=1;rr<=(J+1)/2;rr++)
01659 {
01660 a=nfft_int_2_pow(J-rr);
01661 b=nfft_int_2_pow(rr);
01662
01663 r=NFFT_MIN(rr,J-rr);
01664
01665 n[0]=ths->sigma*nfft_int_2_pow(r);
01666 if(a<b)
01667 n[1]=ths->sigma*nfft_int_2_pow(J-r);
01668 else
01669 n[1]=ths->sigma*nfft_int_2_pow(r);
01670 n[2]=ths->sigma*nfft_int_2_pow(J-r);
01671
01672 ths->set_fftw_plan1[rr] =
01673 fftw_plan_dft(3, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
01674 FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
01675 ths->set_fftw_plan2[rr] =
01676 fftw_plan_dft(3, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
01677 FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
01678 }
01679
01680
01681 for(r=0;r<=nfft_ld(m);r++)
01682 {
01683 N[0]=nfft_int_2_pow(J-r); n[0]=ths->sigma*N[0];
01684 N[1]=nfft_int_2_pow(J-r); n[1]=ths->sigma*N[1];
01685
01686 if(N[0]>m)
01687 {
01688 nfft_init_guru(&(ths->set_nfft_plan_1d[r]),1,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
01689 ths->set_nfft_plan_1d[r].nfft_flags = ths->set_nfft_plan_1d[r].nfft_flags | FG_PSI;
01690 ths->set_nfft_plan_1d[r].K=ths->act_nfft_plan->K;
01691 ths->set_nfft_plan_1d[r].psi=ths->act_nfft_plan->psi;
01692 nfft_init_guru(&(ths->set_nfft_plan_2d[r]),2,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
01693 ths->set_nfft_plan_2d[r].nfft_flags = ths->set_nfft_plan_2d[r].nfft_flags | FG_PSI;
01694 ths->set_nfft_plan_2d[r].K=ths->act_nfft_plan->K;
01695 ths->set_nfft_plan_2d[r].psi=ths->act_nfft_plan->psi;
01696 }
01697 }
01698
01699
01700
01701 N[0]=nfft_int_2_pow(J/2+1); n[0]=ths->sigma*N[0];
01702 N[1]=nfft_int_2_pow(J/2+1); n[1]=ths->sigma*N[1];
01703 N[2]=nfft_int_2_pow(J/2+1); n[2]=ths->sigma*N[2];
01704 nfft_init_guru(ths->center_nfft_plan,3,N,M,n, m, MALLOC_F| FFTW_INIT,
01705 FFTW_MEASURE);
01706 ths->center_nfft_plan->x= ths->act_nfft_plan->x;
01707 ths->center_nfft_plan->nfft_flags = ths->center_nfft_plan->nfft_flags|
01708 FG_PSI;
01709 ths->center_nfft_plan->K=ths->act_nfft_plan->K;
01710 ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
01711
01712 if(ths->flags & NSDFT)
01713 {
01714 ths->index_sparse_to_full=(int*)fftw_malloc(ths->N_total*sizeof(int));
01715 init_index_sparse_to_full_3d(ths);
01716 }
01717 }
01718
01719 #ifdef GAUSSIAN
01720 void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
01721 {
01722 ths->d=d;
01723
01724 if(ths->d==2)
01725 nsfft_init_2d(ths, J, M, m, flags);
01726 else
01727 nsfft_init_3d(ths, J, M, m, flags);
01728 }
01729 #else
01730 void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
01731 {
01732 fprintf(stderr,
01733 "\nError in kernel/nsfft_init: require GAUSSIAN window function\n");
01734 }
01735 #endif
01736
01737 void nsfft_finalize_2d(nsfft_plan *ths)
01738 {
01739 int r;
01740
01741 if(ths->flags & NSDFT)
01742 fftw_free(ths->index_sparse_to_full);
01743
01744
01745 ths->center_nfft_plan->nfft_flags = ths->center_nfft_plan->nfft_flags ^ FG_PSI;
01746 nfft_finalize(ths->center_nfft_plan);
01747
01748
01749 for(r=0;r<=nfft_ld(ths->act_nfft_plan->m);r++)
01750 {
01751 ths->set_nfft_plan_1d[r].nfft_flags =
01752 ths->set_nfft_plan_1d[r].nfft_flags ^ FG_PSI;
01753 nfft_finalize(&(ths->set_nfft_plan_1d[r]));
01754 }
01755
01756
01757 ths->act_nfft_plan->my_fftw_plan2=ths->set_fftw_plan2[0];
01758 ths->act_nfft_plan->my_fftw_plan1=ths->set_fftw_plan1[0];
01759
01760 for(r=1;r<=ths->J/2;r++)
01761 {
01762 fftw_destroy_plan(ths->set_fftw_plan2[r]);
01763 fftw_destroy_plan(ths->set_fftw_plan1[r]);
01764 }
01765
01766
01767 nfft_finalize(ths->act_nfft_plan);
01768
01769 fftw_free(ths->set_nfft_plan_1d);
01770
01771 fftw_free(ths->set_fftw_plan2);
01772 fftw_free(ths->set_fftw_plan1);
01773
01774 fftw_free(ths->x_transposed);
01775
01776 fftw_free(ths->f_hat);
01777 fftw_free(ths->f);
01778 }
01779
01780 void nsfft_finalize_3d(nsfft_plan *ths)
01781 {
01782 int r;
01783
01784 if(ths->flags & NSDFT)
01785 fftw_free(ths->index_sparse_to_full);
01786
01787
01788 ths->center_nfft_plan->nfft_flags = ths->center_nfft_plan->nfft_flags ^ FG_PSI;
01789 nfft_finalize(ths->center_nfft_plan);
01790
01791
01792 for(r=0;r<=nfft_ld(ths->act_nfft_plan->m);r++)
01793 {
01794 if(nfft_int_2_pow(ths->J-r)>ths->act_nfft_plan->m)
01795 {
01796 ths->set_nfft_plan_2d[r].nfft_flags = ths->set_nfft_plan_2d[r].nfft_flags ^ FG_PSI;
01797 nfft_finalize(&(ths->set_nfft_plan_2d[r]));
01798 ths->set_nfft_plan_1d[r].nfft_flags = ths->set_nfft_plan_1d[r].nfft_flags ^ FG_PSI;
01799 nfft_finalize(&(ths->set_nfft_plan_1d[r]));
01800 }
01801 }
01802
01803
01804 ths->act_nfft_plan->my_fftw_plan2=ths->set_fftw_plan2[0];
01805 ths->act_nfft_plan->my_fftw_plan1=ths->set_fftw_plan1[0];
01806
01807 for(r=1;r<=(ths->J+1)/2;r++)
01808 {
01809 fftw_destroy_plan(ths->set_fftw_plan2[r]);
01810 fftw_destroy_plan(ths->set_fftw_plan1[r]);
01811 }
01812
01813
01814 nfft_finalize(ths->act_nfft_plan);
01815
01816 fftw_free(ths->set_nfft_plan_1d);
01817 fftw_free(ths->set_nfft_plan_2d);
01818
01819 fftw_free(ths->set_fftw_plan2);
01820 fftw_free(ths->set_fftw_plan1);
01821
01822 fftw_free(ths->x_102);
01823 fftw_free(ths->x_201);
01824 fftw_free(ths->x_120);
01825 fftw_free(ths->x_021);
01826
01827 fftw_free(ths->f_hat);
01828 fftw_free(ths->f);
01829 }
01830
01831 void nsfft_finalize(nsfft_plan *ths)
01832 {
01833 if(ths->d==2)
01834 nsfft_finalize_2d(ths);
01835 else
01836 nsfft_finalize_3d(ths);
01837 }