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