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

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409