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