00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "config.h"
00022
00023 #include <stdio.h>
00024 #include <math.h>
00025 #include <string.h>
00026 #include <stdlib.h>
00027 #ifdef HAVE_COMPLEX_H
00028 #include <complex.h>
00029 #endif
00030 #include "nfft3util.h"
00031 #include "nfft3.h"
00032 #include "infft.h"
00033
00034
00035 #define MACRO_nndft_init_result_trafo memset(f,0,ths->M_total*sizeof(double _Complex));
00036 #define MACRO_nndft_init_result_conjugated MACRO_nndft_init_result_trafo
00037 #define MACRO_nndft_init_result_adjoint memset(f_hat,0,ths->N_total*sizeof(double _Complex));
00038 #define MACRO_nndft_init_result_transposed MACRO_nndft_init_result_adjoint
00039
00040 #define MACRO_nndft_sign_trafo (-2.0*PI)
00041 #define MACRO_nndft_sign_conjugated (+2.0*PI)
00042 #define MACRO_nndft_sign_adjoint (+2.0*PI)
00043 #define MACRO_nndft_sign_transposed (-2.0*PI)
00044
00045 #define MACRO_nndft_compute_trafo (*fj) += (*f_hat_k)*cexp(+ _Complex_I*omega);
00046
00047 #define MACRO_nndft_compute_conjugated MACRO_nndft_compute_trafo
00048
00049 #define MACRO_nndft_compute_adjoint (*f_hat_k) += (*fj)*cexp(+ _Complex_I*omega);
00050
00051 #define MACRO_nndft_compute_transposed MACRO_nndft_compute_adjoint
00052
00053 #define MACRO_nndft(which_one) \
00054 void nnfft_ ## which_one ## _direct (nnfft_plan *ths) \
00055 { \
00056 int j; \
00057 int t; \
00058 int l; \
00059 double _Complex *f_hat, *f; \
00060 double _Complex *f_hat_k; \
00061 double _Complex *fj; \
00062 double omega; \
00063 \
00064 f_hat=ths->f_hat; f=ths->f; \
00065 \
00066 MACRO_nndft_init_result_ ## which_one \
00067 \
00068 for(j=0, fj=f; j<ths->M_total; j++, fj++) \
00069 { \
00070 for(l=0, f_hat_k=f_hat; l<ths->N_total; l++, f_hat_k++) \
00071 { \
00072 omega=0.0; \
00073 for(t = 0; t<ths->d; t++) \
00074 omega+=ths->v[l*ths->d+t] * ths->x[j*ths->d+t] * ths->N[t]; \
00075 \
00076 omega*= MACRO_nndft_sign_ ## which_one; \
00077 \
00078 MACRO_nndft_compute_ ## which_one \
00079 \
00080 } \
00081 } \
00082 } \
00083
00084 MACRO_nndft(trafo)
00085 MACRO_nndft(adjoint)
00086
00089 static void nnfft_uo(nnfft_plan *ths,int j,int *up,int *op,int act_dim)
00090 {
00091 double c;
00092 int u,o;
00093
00094 c = ths->v[j*ths->d+act_dim] * ths->n[act_dim];
00095
00096 u = c; o = c;
00097 if(c < 0)
00098 u = u-1;
00099 else
00100 o = o+1;
00101
00102 u = u - (ths->m); o = o + (ths->m);
00103
00104 up[0]=u; op[0]=o;
00105 }
00106
00110 #define MACRO_nnfft_B_init_result_A memset(f,0,ths->N_total*sizeof(double _Complex));
00111 #define MACRO_nnfft_B_init_result_T memset(g,0,ths->aN1_total*sizeof(double _Complex));
00112
00113 #define MACRO_nnfft_B_PRE_FULL_PSI_compute_A { \
00114 (*fj) += ths->psi[ix] * g[ths->psi_index_g[ix]]; \
00115 }
00116
00117 #define MACRO_nnfft_B_PRE_FULL_PSI_compute_T { \
00118 g[ths->psi_index_g[ix]] += ths->psi[ix] * (*fj); \
00119 }
00120
00121 #define MACRO_nnfft_B_compute_A { \
00122 (*fj) += phi_prod[ths->d] * g[ll_plain[ths->d]]; \
00123 }
00124
00125 #define MACRO_nnfft_B_compute_T { \
00126 g[ll_plain[ths->d]] += phi_prod[ths->d] * (*fj); \
00127 }
00128
00129 #define MACRO_with_PRE_LIN_PSI (ths->psi[(ths->K+1)*t2+y_u[t2]]* \
00130 (y_u[t2]+1-y[t2]) + \
00131 ths->psi[(ths->K+1)*t2+y_u[t2]+1]* \
00132 (y[t2]-y_u[t2]))
00133 #define MACRO_with_PRE_PSI ths->psi[(j*ths->d+t2)*(2*ths->m+2)+lj[t2]]
00134 #define MACRO_without_PRE_PSI PHI(-ths->v[j*ths->d+t2]+ \
00135 ((double)l[t2])/ths->N1[t2], t2)
00136
00137 #define MACRO_init_uo_l_lj_t { \
00138 for(t = ths->d-1; t>=0; t--) \
00139 { \
00140 nnfft_uo(ths,j,&u[t],&o[t],t); \
00141 l[t] = u[t]; \
00142 lj[t] = 0; \
00143 } \
00144 t++; \
00145 }
00146
00147 #define MACRO_update_with_PRE_PSI_LIN { \
00148 for(t2=t; t2<ths->d; t2++) \
00149 { \
00150 y[t2] = fabs(((-ths->N1[t2]*ths->v[j*ths->d+t2]+(double)l[t2]) \
00151 * ((double)ths->K))/(ths->m+1)); \
00152 y_u[t2] = (int)y[t2]; \
00153 } \
00154 }
00155
00156 #define MACRO_update_phi_prod_ll_plain(which_one) { \
00157 for(t2=t; t2<ths->d; t2++) \
00158 { \
00159 phi_prod[t2+1]=phi_prod[t2]* MACRO_ ## which_one; \
00160 ll_plain[t2+1]=ll_plain[t2]*ths->aN1[t2] + \
00161 (l[t2]+ths->aN1[t2]*3/2)%ths->aN1[t2]; \
00162 \
00163 } \
00164 }
00165
00166 #define MACRO_count_uo_l_lj_t { \
00167 for(t = ths->d-1; (t>0)&&(l[t]==o[t]); t--) \
00168 { \
00169 l[t] = u[t]; \
00170 lj[t] = 0; \
00171 } \
00172 \
00173 l[t]++; \
00174 lj[t]++; \
00175 }
00176
00177 #define MACRO_nnfft_B(which_one) \
00178 static inline void nnfft_B_ ## which_one (nnfft_plan *ths) \
00179 { \
00180 int lprod; \
00181 int u[ths->d], o[ths->d]; \
00182 int t, t2; \
00183 int j; \
00184 int l_L, ix; \
00185 int l[ths->d]; \
00186 int lj[ths->d]; \
00187 int ll_plain[ths->d+1]; \
00188 double phi_prod[ths->d+1]; \
00189 double _Complex *f, *g; \
00190 double _Complex *fj; \
00191 double y[ths->d]; \
00192 int y_u[ths->d]; \
00193 \
00194 f=ths->f_hat; g=ths->F; \
00195 \
00196 MACRO_nnfft_B_init_result_ ## which_one \
00197 \
00198 if(ths->nnfft_flags & PRE_FULL_PSI) \
00199 { \
00200 for(ix=0, j=0, fj=f; j<ths->N_total; j++,fj++) \
00201 for(l_L=0; l_L<ths->psi_index_f[j]; l_L++, ix++) \
00202 MACRO_nnfft_B_PRE_FULL_PSI_compute_ ## which_one; \
00203 return; \
00204 } \
00205 \
00206 phi_prod[0]=1; \
00207 ll_plain[0]=0; \
00208 \
00209 for(t=0,lprod = 1; t<ths->d; t++) \
00210 lprod *= (2*ths->m+2); \
00211 \
00212 if(ths->nnfft_flags & PRE_PSI) \
00213 { \
00214 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
00215 { \
00216 MACRO_init_uo_l_lj_t; \
00217 \
00218 for(l_L=0; l_L<lprod; l_L++) \
00219 { \
00220 MACRO_update_phi_prod_ll_plain(with_PRE_PSI); \
00221 \
00222 MACRO_nnfft_B_compute_ ## which_one; \
00223 \
00224 MACRO_count_uo_l_lj_t; \
00225 } \
00226 } \
00227 return; \
00228 } \
00229 \
00230 if(ths->nnfft_flags & PRE_LIN_PSI) \
00231 { \
00232 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
00233 { \
00234 MACRO_init_uo_l_lj_t; \
00235 \
00236 for(l_L=0; l_L<lprod; l_L++) \
00237 { \
00238 MACRO_update_with_PRE_PSI_LIN; \
00239 \
00240 MACRO_update_phi_prod_ll_plain(with_PRE_LIN_PSI); \
00241 \
00242 MACRO_nnfft_B_compute_ ## which_one; \
00243 \
00244 MACRO_count_uo_l_lj_t; \
00245 } \
00246 } \
00247 return; \
00248 } \
00249 \
00250 \
00251 for(j=0, fj=f; j<ths->N_total; j++, fj++) \
00252 { \
00253 \
00254 MACRO_init_uo_l_lj_t; \
00255 \
00256 for(l_L=0; l_L<lprod; l_L++) \
00257 { \
00258 MACRO_update_phi_prod_ll_plain(without_PRE_PSI); \
00259 \
00260 MACRO_nnfft_B_compute_ ## which_one; \
00261 \
00262 MACRO_count_uo_l_lj_t; \
00263 } \
00264 } \
00265 }
00266
00267 MACRO_nnfft_B(A)
00268 MACRO_nnfft_B(T)
00269
00270 static inline void nnfft_D (nnfft_plan *ths){
00271 int j,t;
00272 double tmp;
00273
00274 if(ths->nnfft_flags & PRE_PHI_HUT)
00275 {
00276 for(j=0; j<ths->M_total; j++)
00277 ths->f[j] *= ths->c_phi_inv[j];
00278 }
00279 else
00280 {
00281 for(j=0; j<ths->M_total; j++)
00282 {
00283 tmp = 1.0;
00284
00285 for(t=0; t<ths->d; t++)
00286 tmp*= 1.0 /((PHI_HUT(ths->x[ths->d*j + t]*((double)ths->N[t]),t)) );
00287 ths->f[j] *= tmp;
00288 }
00289 }
00290 }
00291
00294 void nnfft_trafo(nnfft_plan *ths)
00295 {
00296 int j,t;
00297
00298 nnfft_B_T(ths);
00299
00300 for(j=0;j<ths->M_total;j++) {
00301 for(t=0;t<ths->d;t++) {
00302 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] / ((double)ths->sigma[t]);
00303 }
00304 }
00305
00306
00307 ths->direct_plan->f = ths->f;
00308
00309 nfft_trafo(ths->direct_plan);
00310
00311 for(j=0;j<ths->M_total;j++) {
00312 for(t=0;t<ths->d;t++) {
00313 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] * ((double)ths->sigma[t]);
00314 }
00315 }
00316
00317 nnfft_D(ths);
00318 }
00319
00320 void nnfft_adjoint(nnfft_plan *ths)
00321 {
00322 int j,t;
00323
00324 nnfft_D(ths);
00325
00326 for(j=0;j<ths->M_total;j++) {
00327 for(t=0;t<ths->d;t++) {
00328 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] / ((double)ths->sigma[t]);
00329 }
00330 }
00331
00332
00333 ths->direct_plan->f=ths->f;
00334
00335 nfft_adjoint(ths->direct_plan);
00336
00337 for(j=0;j<ths->M_total;j++) {
00338 for(t=0;t<ths->d;t++) {
00339 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] * ((double)ths->sigma[t]);
00340 }
00341 }
00342
00343 nnfft_B_A(ths);
00344 }
00345
00348 void nnfft_precompute_phi_hut(nnfft_plan *ths)
00349 {
00350 int j;
00351 int t;
00352 double tmp;
00353
00354 ths->c_phi_inv= (double*)nfft_malloc(ths->M_total*sizeof(double));
00355
00356 for(j=0; j<ths->M_total; j++)
00357 {
00358 tmp = 1.0;
00359 for(t=0; t<ths->d; t++)
00360 tmp*= 1.0 /(PHI_HUT(ths->x[ths->d*j + t]*((double)ths->N[t]),t));
00361 ths->c_phi_inv[j]=tmp;
00362 }
00363 }
00364
00365
00368 void nnfft_precompute_lin_psi(nnfft_plan *ths)
00369 {
00370 int t;
00371 int j;
00372 double step;
00374 nfft_precompute_lin_psi(ths->direct_plan);
00375
00376 for (t=0; t<ths->d; t++)
00377 {
00378 step=((double)(ths->m+1))/(ths->K*ths->N1[t]);
00379 for(j=0;j<=ths->K;j++)
00380 {
00381 ths->psi[(ths->K+1)*t + j] = PHI(j*step,t);
00382 }
00383 }
00384 }
00385
00386 void nnfft_precompute_psi(nnfft_plan *ths)
00387 {
00388 int t;
00389 int j;
00390 int l;
00391 int lj;
00392 int u, o;
00394 for (t=0; t<ths->d; t++)
00395 for(j=0;j<ths->N_total;j++)
00396 {
00397 nnfft_uo(ths,j,&u,&o,t);
00398
00399 for(l=u, lj=0; l <= o; l++, lj++)
00400 ths->psi[(j*ths->d+t)*(2*ths->m+2)+lj]=
00401 (PHI((-ths->v[j*ths->d+t]+((double)l)/((double)ths->N1[t])),t));
00402 }
00403
00404 for(j=0;j<ths->M_total;j++) {
00405 for(t=0;t<ths->d;t++) {
00406 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] / ((double)ths->sigma[t]);
00407 }
00408 }
00409
00410 nfft_precompute_psi(ths->direct_plan);
00411
00412 for(j=0;j<ths->M_total;j++) {
00413 for(t=0;t<ths->d;t++) {
00414 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] * ((double)ths->sigma[t]);
00415 }
00416 }
00417
00418 }
00419
00420
00421
00425 void nnfft_precompute_full_psi(nnfft_plan *ths)
00426 {
00427 int t,t2;
00428 int j;
00429 int l_L;
00430 int l[ths->d];
00431 int lj[ths->d];
00432 int ll_plain[ths->d+1];
00433 int lprod;
00434 int u[ths->d], o[ths->d];
00436 double phi_prod[ths->d+1];
00437
00438 int ix,ix_old;
00439
00440 for(j=0;j<ths->M_total;j++) {
00441 for(t=0;t<ths->d;t++) {
00442 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] / ((double)ths->sigma[t]);
00443 }
00444 }
00445
00446 nnfft_precompute_psi(ths);
00447
00448 nfft_precompute_full_psi(ths->direct_plan);
00449
00450 for(j=0;j<ths->M_total;j++) {
00451 for(t=0;t<ths->d;t++) {
00452 ths->x[j*ths->d+t]= ths->x[j*ths->d+t] * ((double)ths->sigma[t]);
00453 }
00454 }
00455
00456 phi_prod[0]=1;
00457 ll_plain[0]=0;
00458
00459 for(t=0,lprod = 1; t<ths->d; t++)
00460 lprod *= 2*ths->m+2;
00461
00462 for(j=0,ix=0,ix_old=0; j<ths->N_total; j++)
00463 {
00464 MACRO_init_uo_l_lj_t;
00465
00466 for(l_L=0; l_L<lprod; l_L++, ix++)
00467 {
00468 MACRO_update_phi_prod_ll_plain(without_PRE_PSI);
00469
00470 ths->psi_index_g[ix]=ll_plain[ths->d];
00471 ths->psi[ix]=phi_prod[ths->d];
00472
00473 MACRO_count_uo_l_lj_t;
00474 }
00475
00476
00477 ths->psi_index_f[j]=ix-ix_old;
00478 ix_old=ix;
00479 }
00480 }
00481
00482 static void nnfft_init_help(nnfft_plan *ths, int m2, unsigned nfft_flags, unsigned fftw_flags)
00483 {
00484 int t;
00485 int lprod;
00486 int N2[ths->d];
00487
00488 ths->aN1 = (int*) nfft_malloc(ths->d*sizeof(int));
00489
00490 ths->a = (double*) nfft_malloc(ths->d*sizeof(double));
00491
00492 ths->sigma = (double*) nfft_malloc(ths->d*sizeof(double));
00493
00494 ths->n = ths->N1;
00495
00496 ths->aN1_total=1;
00497
00498 for(t = 0; t<ths->d; t++) {
00499 ths->a[t] = 1.0 + (2.0*((double)ths->m))/((double)ths->N1[t]);
00500 ths->aN1[t] = ths->a[t] * ((double)ths->N1[t]);
00501
00502 if(ths->aN1[t]%2 != 0)
00503 ths->aN1[t] = ths->aN1[t] +1;
00504
00505 ths->aN1_total*=ths->aN1[t];
00506 ths->sigma[t] = ((double) ths->N1[t] )/((double) ths->N[t]);;
00507
00508
00509 N2[t] = ceil(ths->sigma[t]*(ths->aN1[t]));
00510
00511
00512 if(N2[t]%2 != 0)
00513 N2[t] = N2[t] +1;
00514 }
00515
00516 WINDOW_HELP_INIT
00517
00518 if(ths->nnfft_flags & MALLOC_X)
00519 ths->x = (double*)nfft_malloc(ths->d*ths->M_total*sizeof(double));
00520 if(ths->nnfft_flags & MALLOC_F)
00521 ths->f=(double _Complex*)nfft_malloc(ths->M_total*sizeof(double _Complex));
00522
00523 if(ths->nnfft_flags & MALLOC_V)
00524 ths->v = (double*)nfft_malloc(ths->d*ths->N_total*sizeof(double));
00525 if(ths->nnfft_flags & MALLOC_F_HAT)
00526 ths->f_hat = (double _Complex*)nfft_malloc(ths->N_total*sizeof(double _Complex));
00527
00528 if(ths->nnfft_flags & PRE_LIN_PSI)
00529 {
00530 ths->K=(1U<< 10)*(ths->m+1);
00531 ths->psi = (double*) nfft_malloc((ths->K+1)*ths->d*sizeof(double));
00532 }
00533
00534 if(ths->nnfft_flags & PRE_PSI)
00535 ths->psi = (double*)nfft_malloc(ths->N_total*ths->d*(2*ths->m+2)*sizeof(double));
00536
00537 if(ths->nnfft_flags & PRE_FULL_PSI)
00538 {
00539 for(t=0,lprod = 1; t<ths->d; t++)
00540 lprod *= 2*ths->m+2;
00541
00542 ths->psi = (double*)nfft_malloc(ths->N_total*lprod*sizeof(double));
00543
00544 ths->psi_index_f = (int*) nfft_malloc(ths->N_total*sizeof(int));
00545 ths->psi_index_g = (int*) nfft_malloc(ths->N_total*lprod*sizeof(int));
00546 }
00547
00548 ths->direct_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
00549
00550 nfft_init_guru(ths->direct_plan, ths->d, ths->aN1, ths->M_total, N2, m2,
00551 nfft_flags, fftw_flags);
00552
00553 ths->direct_plan->x = ths->x;
00554 ths->direct_plan->f = ths->f;
00555 ths->F = ths->direct_plan->f_hat;
00556
00557 ths->mv_trafo = (void (*) (void* ))nnfft_trafo;
00558 ths->mv_adjoint = (void (*) (void* ))nnfft_adjoint;
00559 }
00560
00561 void nnfft_init_guru(nnfft_plan *ths, int d, int N_total, int M_total, int *N, int *N1,
00562 int m, unsigned nnfft_flags)
00563 {
00564 int t;
00566 unsigned nfft_flags;
00567 unsigned fftw_flags;
00568
00569 ths->d= d;
00570 ths->M_total= M_total;
00571 ths->N_total= N_total;
00572 ths->m= m;
00573 ths->nnfft_flags= nnfft_flags;
00574 fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
00575 nfft_flags= PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE;
00576
00577 if(ths->nnfft_flags & PRE_PSI)
00578 nfft_flags = nfft_flags | PRE_PSI;
00579
00580 if(ths->nnfft_flags & PRE_FULL_PSI)
00581 nfft_flags = nfft_flags | PRE_FULL_PSI;
00582
00583 if(ths->nnfft_flags & PRE_LIN_PSI)
00584 nfft_flags = nfft_flags | PRE_LIN_PSI;
00585
00586 ths->N = (int*) nfft_malloc(ths->d*sizeof(int));
00587 ths->N1 = (int*) nfft_malloc(ths->d*sizeof(int));
00588
00589 for(t=0; t<d; t++) {
00590 ths->N[t] = N[t];
00591 ths->N1[t] = N1[t];
00592 }
00593 nnfft_init_help(ths,m,nfft_flags,fftw_flags);
00594 }
00595
00596 void nnfft_init(nnfft_plan *ths, int d, int N_total, int M_total, int *N)
00597 {
00598 int t;
00600 unsigned nfft_flags;
00601 unsigned fftw_flags;
00602
00603 ths->d = d;
00604 ths->M_total = M_total;
00605 ths->N_total = N_total;
00606
00607
00608
00609
00610
00611
00612
00613 ths->N = (int*) nfft_malloc(ths->d*sizeof(int));
00614 ths->N1 = (int*) nfft_malloc(ths->d*sizeof(int));
00615
00616 for(t=0; t<d; t++) {
00617 ths->N[t] = N[t];
00618
00619
00620 ths->N1[t] = ceil(1.5*ths->N[t]);
00621
00622
00623 if(ths->N1[t]%2 != 0)
00624 ths->N1[t] = ths->N1[t] +1;
00625 }
00626
00627 ths->nnfft_flags=PRE_PSI| PRE_PHI_HUT| MALLOC_X| MALLOC_V| MALLOC_F_HAT| MALLOC_F;
00628 nfft_flags= PRE_PSI| PRE_PHI_HUT| MALLOC_F_HAT| FFTW_INIT| FFT_OUT_OF_PLACE;
00629
00630 fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;
00631
00632 nnfft_init_help(ths,ths->m,nfft_flags,fftw_flags);
00633 }
00634
00635 void nnfft_finalize(nnfft_plan *ths)
00636 {
00637
00638 nfft_finalize(ths->direct_plan);
00639
00640 nfft_free(ths->direct_plan);
00641
00642 nfft_free(ths->aN1);
00643 nfft_free(ths->N);
00644 nfft_free(ths->N1);
00645
00646 if(ths->nnfft_flags & PRE_FULL_PSI)
00647 {
00648 nfft_free(ths->psi_index_g);
00649 nfft_free(ths->psi_index_f);
00650 nfft_free(ths->psi);
00651 }
00652
00653 if(ths->nnfft_flags & PRE_PSI)
00654 nfft_free(ths->psi);
00655
00656 if(ths->nnfft_flags & PRE_LIN_PSI)
00657 nfft_free(ths->psi);
00658
00659 if(ths->nnfft_flags & PRE_PHI_HUT)
00660 nfft_free(ths->c_phi_inv);
00661
00662 if(ths->nnfft_flags & MALLOC_F)
00663 nfft_free(ths->f);
00664
00665 if(ths->nnfft_flags & MALLOC_F_HAT)
00666 nfft_free(ths->f_hat);
00667
00668 if(ths->nnfft_flags & MALLOC_X)
00669 nfft_free(ths->x);
00670
00671 if(ths->nnfft_flags & MALLOC_V)
00672 nfft_free(ths->v);
00673 }