00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00029 #include <math.h>
00030 #include <stdlib.h>
00031 #include <complex.h>
00032
00033 #include "nfft3util.h"
00034 #include "nfft3.h"
00035
00042 double GLOBAL_elapsed_time;
00043
00047 static int linogram_grid(int T, int R, double *x, double *w)
00048 {
00049 int t, r;
00050 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00051
00052 for(t=-T/2; t<T/2; t++)
00053 {
00054 for(r=-R/2; r<R/2; r++)
00055 {
00056 if(t<0)
00057 {
00058 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
00059 x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
00060 }
00061 else
00062 {
00063 x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
00064 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
00065 }
00066 if (r==0)
00067 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00068 else
00069 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00070 }
00071 }
00072
00073 return T*R;
00074 }
00075
00077 static int linogram_dft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00078 {
00079 int j,k;
00080 nfft_plan my_nfft_plan;
00082 double *x, *w;
00084 int N[2],n[2];
00085 int M=T*R;
00087 N[0]=NN; n[0]=2*N[0];
00088 N[1]=NN; n[1]=2*N[1];
00090 x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00091 if (x==NULL)
00092 return -1;
00093
00094 w = (double *)nfft_malloc(T*R*(sizeof(double)));
00095 if (w==NULL)
00096 return -1;
00097
00099 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00100 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00101 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00102
00104 linogram_grid(T,R,x,w);
00105 for(j=0;j<my_nfft_plan.M_total;j++)
00106 {
00107 my_nfft_plan.x[2*j+0] = x[2*j+0];
00108 my_nfft_plan.x[2*j+1] = x[2*j+1];
00109 }
00110
00111
00113 for(k=0;k<my_nfft_plan.N_total;k++)
00114 my_nfft_plan.f_hat[k] = f_hat[k];
00115
00117 GLOBAL_elapsed_time=nfft_second();
00118 ndft_trafo(&my_nfft_plan);
00119 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00120
00122 for(j=0;j<my_nfft_plan.M_total;j++)
00123 f[j] = my_nfft_plan.f[j];
00124
00126 nfft_finalize(&my_nfft_plan);
00127 nfft_free(x);
00128 nfft_free(w);
00129
00130 return EXIT_SUCCESS;
00131 }
00132
00134 static int linogram_fft(fftw_complex *f_hat, int NN, fftw_complex *f, int T, int R, int m)
00135 {
00136 int j,k;
00137 nfft_plan my_nfft_plan;
00139 double *x, *w;
00141 int N[2],n[2];
00142 int M=T*R;
00144 N[0]=NN; n[0]=2*N[0];
00145 N[1]=NN; n[1]=2*N[1];
00147 x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00148 if (x==NULL)
00149 return -1;
00150
00151 w = (double *)nfft_malloc(T*R*(sizeof(double)));
00152 if (w==NULL)
00153 return -1;
00154
00156 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00157 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00158 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00159
00161 linogram_grid(T,R,x,w);
00162 for(j=0;j<my_nfft_plan.M_total;j++)
00163 {
00164 my_nfft_plan.x[2*j+0] = x[2*j+0];
00165 my_nfft_plan.x[2*j+1] = x[2*j+1];
00166 }
00167
00169 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00170 nfft_precompute_lin_psi(&my_nfft_plan);
00171
00172 if(my_nfft_plan.nfft_flags & PRE_PSI)
00173 nfft_precompute_psi(&my_nfft_plan);
00174
00175 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00176 nfft_precompute_full_psi(&my_nfft_plan);
00177
00179 for(k=0;k<my_nfft_plan.N_total;k++)
00180 my_nfft_plan.f_hat[k] = f_hat[k];
00181
00183 GLOBAL_elapsed_time=nfft_second();
00184 nfft_trafo(&my_nfft_plan);
00185 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00186
00188 for(j=0;j<my_nfft_plan.M_total;j++)
00189 f[j] = my_nfft_plan.f[j];
00190
00192 nfft_finalize(&my_nfft_plan);
00193 nfft_free(x);
00194 nfft_free(w);
00195
00196 return EXIT_SUCCESS;
00197 }
00198
00200 static int inverse_linogram_fft(fftw_complex *f, int T, int R, fftw_complex *f_hat, int NN, int max_i, int m)
00201 {
00202 int j,k;
00203 nfft_plan my_nfft_plan;
00204 solver_plan_complex my_infft_plan;
00206 double *x, *w;
00207 int l;
00209 int N[2],n[2];
00210 int M=T*R;
00212 N[0]=NN; n[0]=2*N[0];
00213 N[1]=NN; n[1]=2*N[1];
00215 x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00216 if (x==NULL)
00217 return -1;
00218
00219 w = (double *)nfft_malloc(T*R*(sizeof(double)));
00220 if (w==NULL)
00221 return -1;
00222
00224 nfft_init_guru(&my_nfft_plan, 2, N, M, n, m,
00225 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00226 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00227
00229 solver_init_advanced_complex(&my_infft_plan,(mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT );
00230
00232 linogram_grid(T,R,x,w);
00233 for(j=0;j<my_nfft_plan.M_total;j++)
00234 {
00235 my_nfft_plan.x[2*j+0] = x[2*j+0];
00236 my_nfft_plan.x[2*j+1] = x[2*j+1];
00237 my_infft_plan.y[j] = f[j];
00238 my_infft_plan.w[j] = w[j];
00239 }
00240
00242 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00243 nfft_precompute_lin_psi(&my_nfft_plan);
00244
00245 if(my_nfft_plan.nfft_flags & PRE_PSI)
00246 nfft_precompute_psi(&my_nfft_plan);
00247
00248 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00249 nfft_precompute_full_psi(&my_nfft_plan);
00250
00252 if(my_infft_plan.flags & PRECOMPUTE_DAMP)
00253 for(j=0;j<my_nfft_plan.N[0];j++)
00254 for(k=0;k<my_nfft_plan.N[1];k++)
00255 {
00256 my_infft_plan.w_hat[j*my_nfft_plan.N[1]+k]=
00257 (sqrt(pow(j-my_nfft_plan.N[0]/2,2)+pow(k-my_nfft_plan.N[1]/2,2))>(my_nfft_plan.N[0]/2)?0:1);
00258 }
00259
00261 for(k=0;k<my_nfft_plan.N_total;k++)
00262 my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
00263
00264 GLOBAL_elapsed_time=nfft_second();
00266 solver_before_loop_complex(&my_infft_plan);
00267
00268 if (max_i<1)
00269 {
00270 l=1;
00271 for(k=0;k<my_nfft_plan.N_total;k++)
00272 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00273 }
00274 else
00275 {
00276 for(l=1;l<=max_i;l++)
00277 {
00278 solver_loop_one_step_complex(&my_infft_plan);
00279 }
00280 }
00281
00282 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00283
00285 for(k=0;k<my_nfft_plan.N_total;k++)
00286 f_hat[k] = my_infft_plan.f_hat_iter[k];
00287
00289 solver_finalize_complex(&my_infft_plan);
00290 nfft_finalize(&my_nfft_plan);
00291 nfft_free(x);
00292 nfft_free(w);
00293
00294 return EXIT_SUCCESS;
00295 }
00296
00298 int comparison_fft(FILE *fp, int N, int T, int R)
00299 {
00300 fftw_plan my_fftw_plan;
00301 fftw_complex *f_hat,*f;
00302 int m,k;
00303 double t_fft, t_dft_linogram;
00304
00305 f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00306 f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*(T*R/4)*5);
00307
00308 my_fftw_plan = fftw_plan_dft_2d(N,N,f_hat,f,FFTW_BACKWARD,FFTW_MEASURE);
00309
00310 for(k=0; k<N*N; k++)
00311 f_hat[k] = (((double)rand())/RAND_MAX) + _Complex_I* (((double)rand())/RAND_MAX);
00312
00313 GLOBAL_elapsed_time=nfft_second();
00314 for(m=0;m<65536/N;m++)
00315 {
00316 fftw_execute(my_fftw_plan);
00317
00318 f_hat[2]=2*f_hat[0];
00319 }
00320 GLOBAL_elapsed_time=nfft_second()-GLOBAL_elapsed_time;
00321 t_fft=N*GLOBAL_elapsed_time/65536;
00322
00323 if(N<256)
00324 {
00325 linogram_dft(f_hat,N,f,T,R,m);
00326 t_dft_linogram=GLOBAL_elapsed_time;
00327 }
00328
00329 for (m=3; m<=9; m+=3)
00330 {
00331 if((m==3)&&(N<256))
00332 fprintf(fp,"%d\t&\t&\t%1.1e&\t%1.1e&\t%d\t",N,t_fft,t_dft_linogram,m);
00333 else
00334 if(m==3)
00335 fprintf(fp,"%d\t&\t&\t%1.1e&\t &\t%d\t",N,t_fft,m);
00336 else
00337 fprintf(fp," \t&\t&\t &\t &\t%d\t",m);
00338
00339 printf("N=%d\tt_fft=%1.1e\tt_dft_linogram=%1.1e\tm=%d\t",N,t_fft,t_dft_linogram,m);
00340
00341 linogram_fft(f_hat,N,f,T,R,m);
00342 fprintf(fp,"%1.1e&\t",GLOBAL_elapsed_time);
00343 printf("t_linogram=%1.1e\t",GLOBAL_elapsed_time);
00344 inverse_linogram_fft(f,T,R,f_hat,N,m+3,m);
00345 if(m==9)
00346 fprintf(fp,"%1.1e\\\\\\hline\n",GLOBAL_elapsed_time);
00347 else
00348 fprintf(fp,"%1.1e\\\\\n",GLOBAL_elapsed_time);
00349 printf("t_ilinogram=%1.1e\n",GLOBAL_elapsed_time);
00350 }
00351
00352 fflush(fp);
00353
00354 nfft_free(f);
00355 nfft_free(f_hat);
00356
00357 return EXIT_SUCCESS;
00358 }
00359
00361 int main(int argc,char **argv)
00362 {
00363 int N;
00364 int T, R;
00365 int M;
00366 double *x, *w;
00367 fftw_complex *f_hat, *f, *f_direct, *f_tilde;
00368 int k;
00369 int max_i;
00370 int m;
00371 double temp1, temp2, E_max=0.0;
00372 FILE *fp1, *fp2;
00373 char filename[30];
00374 int logN;
00375
00376 if( argc!=4 )
00377 {
00378 printf("linogram_fft_test N T R \n");
00379 printf("\n");
00380 printf("N linogram FFT of size NxN \n");
00381 printf("T number of slopes \n");
00382 printf("R number of offsets \n");
00383
00385 printf("\nHence, comparison FFTW, linogram FFT and inverse linogram FFT\n");
00386 fp1=fopen("linogram_comparison_fft.dat","w");
00387 if (fp1==NULL)
00388 return(-1);
00389 for (logN=4; logN<=8; logN++)
00390 comparison_fft(fp1,(1U<< logN), 3*(1U<< logN), 3*(1U<< (logN-1)));
00391 fclose(fp1);
00392
00393 exit(-1);
00394 }
00395
00396 N = atoi(argv[1]);
00397 T = atoi(argv[2]);
00398 R = atoi(argv[3]);
00399 printf("N=%d, linogram grid with T=%d, R=%d => ",N,T,R);
00400
00401 x = (double *)nfft_malloc(5*T*R/2*(sizeof(double)));
00402 w = (double *)nfft_malloc(5*T*R/4*(sizeof(double)));
00403
00404 f_hat = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00405 f = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*5*T*R/4);
00406 f_direct = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*5*T*R/4);
00407 f_tilde = (fftw_complex *)nfft_malloc(sizeof(fftw_complex)*N*N);
00408
00410 M=linogram_grid(T,R,x,w); printf("M=%d.\n",M);
00411
00413 fp1=fopen("input_data_r.dat","r");
00414 fp2=fopen("input_data_i.dat","r");
00415 if ((fp1==NULL) || (fp2==NULL))
00416 return(-1);
00417 for(k=0;k<N*N;k++)
00418 {
00419 fscanf(fp1,"%le ",&temp1);
00420 fscanf(fp2,"%le ",&temp2);
00421 f_hat[k]=temp1+_Complex_I*temp2;
00422 }
00423 fclose(fp1);
00424 fclose(fp2);
00425
00427 linogram_dft(f_hat,N,f_direct,T,R,1);
00428
00429
00431 printf("\nTest of the linogram FFT: \n");
00432 fp1=fopen("linogram_fft_error.dat","w+");
00433 for (m=1; m<=12; m++)
00434 {
00436 linogram_fft(f_hat,N,f,T,R,m);
00437
00439 E_max=nfft_error_l_infty_complex(f_direct,f,M);
00440
00441
00442 printf("m=%2d: E_max = %e\n",m,E_max);
00443 fprintf(fp1,"%e\n",E_max);
00444 }
00445 fclose(fp1);
00446
00448 for (m=3; m<=9; m+=3)
00449 {
00450 printf("\nTest of the inverse linogram FFT for m=%d: \n",m);
00451 sprintf(filename,"linogram_ifft_error%d.dat",m);
00452 fp1=fopen(filename,"w+");
00453 for (max_i=0; max_i<=20; max_i+=2)
00454 {
00456 inverse_linogram_fft(f_direct,T,R,f_tilde,N,max_i,m);
00457
00459 E_max=nfft_error_l_infty_complex(f_hat,f_tilde,N*N);
00460 printf("%3d iterations: E_max = %e\n",max_i,E_max);
00461 fprintf(fp1,"%e\n",E_max);
00462 }
00463 fclose(fp1);
00464 }
00465
00467 nfft_free(x);
00468 nfft_free(w);
00469 nfft_free(f_hat);
00470 nfft_free(f);
00471 nfft_free(f_direct);
00472 nfft_free(f_tilde);
00473
00474 return 0;
00475 }
00476