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