NFFT Logo 3.2.2
nnfft.c
00001 /*
00002  * Copyright (c) 2002, 2012 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 /* $Id: nnfft.c 3896 2012-10-10 12:19:26Z tovo $ */
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      } /* for(l) */                                                           \
00081    } /* for(j) */                                                             \
00082 } /* nndft_trafo */                                                           \
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     } /* for(t) */                                                            \
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     } /* for(t2) */                                                           \
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       /* 3/2 because of the (not needed) fftshift and to be in [0 aN1[t2]]?!*/\
00163     } /* for(t2) */                                                           \
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     } /* for(t) */                                                            \
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             } /* for(l_L) */                                                  \
00226         } /* for(j) */                                                        \
00227       return;                                                                 \
00228     } /* if(PRE_PSI) */                                                       \
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             } /* for(l_L) */                                                  \
00246         } /* for(j) */                                                        \
00247       return;                                                                 \
00248     } /* if(PRE_LIN_PSI) */                                                   \
00249                                                                               \
00250   /* no precomputed psi at all */                                             \
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         } /* for(l_L) */                                                      \
00264     } /* for(j) */                                                            \
00265 } /* nnfft_B */
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     /* multiply with N1, because x was modified */
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   /* allows for external swaps of ths->f */
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 } /* nnfft_trafo */
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   /* allows for external swaps of ths->f */
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 } /* nnfft_adjoint */
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 } /* nnfft_phi_hut */
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         } /* for(j) */
00383     } /* for(t) */
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       } /* for(j) */
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   /* for(t) */
00418 } /* nfft_precompute_psi */
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         } /* for(l_L) */
00475 
00476 
00477       ths->psi_index_f[j]=ix-ix_old;
00478       ix_old=ix;
00479     } /* for(j) */
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     /* aN1 should be even */
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     /* take the same oversampling factor in the inner NFFT */
00509     N2[t] = ceil(ths->sigma[t]*(ths->aN1[t]));
00510     
00511     /* N2 should be even */
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   /* m should be greater to get the same accuracy as the nfft */
00608 /* Was soll dieser Ausdruck machen? Es handelt sich um eine Ganzzahl!
00609 
00610   WINDOW_HELP_ESTIMATE_m;
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     /* the standard oversampling factor in the nnfft is 1.5 */
00620     ths->N1[t] = ceil(1.5*ths->N[t]);
00621     
00622     /* N1 should be even */
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409