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