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