00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "config.h"
00021
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #ifdef HAVE_COMPLEX_H
00025 #include <complex.h>
00026 #endif
00027
00028 #include "nfft3util.h"
00029 #include "nfft3.h"
00030
00040 static void reconstruct(char* filename, int N, int M, int weight)
00041 {
00042 int j;
00043 double weights;
00044 double real,imag;
00045 nfft_plan my_plan;
00046 FILE* fin;
00047 FILE* fweight;
00048 FILE *fout_real;
00049 FILE *fout_imag;
00050 int my_N[2],my_n[2];
00051 int flags = PRE_PHI_HUT| PRE_PSI |MALLOC_X| MALLOC_F_HAT|
00052 MALLOC_F| FFTW_INIT| FFT_OUT_OF_PLACE|
00053 FFTW_MEASURE| FFTW_DESTROY_INPUT;
00054
00055
00056 my_N[0]=N; my_n[0]=ceil(N*1.2);
00057 my_N[1]=N; my_n[1]=ceil(N*1.2);
00058 nfft_init_guru(&my_plan, 2, my_N, M, my_n, 6,flags,
00059 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00060
00061 fin=fopen(filename,"r");
00062
00063 fweight=fopen("weights.dat","r");
00064 for(j=0;j<my_plan.M_total;j++)
00065 {
00066 fscanf(fweight,"%le ",&weights);
00067 fscanf(fin,"%le %le %le %le",&my_plan.x[2*j+0],&my_plan.x[2*j+1],&real,&imag);
00068 my_plan.f[j] = real + _Complex_I*imag;
00069 if (weight)
00070 my_plan.f[j] = my_plan.f[j] * weights;
00071 }
00072 fclose(fweight);
00073
00074
00075 if(my_plan.nfft_flags & PRE_PSI)
00076 nfft_precompute_psi(&my_plan);
00077
00078
00079 if(my_plan.nfft_flags & PRE_FULL_PSI)
00080 nfft_precompute_full_psi(&my_plan);
00081
00082
00083
00084 nfft_adjoint(&my_plan);
00085
00086 fout_real=fopen("output_real.dat","w");
00087 fout_imag=fopen("output_imag.dat","w");
00088
00089 for (j=0;j<N*N;j++) {
00090 fprintf(fout_real,"%le ",creal(my_plan.f_hat[j]));
00091 fprintf(fout_imag,"%le ",cimag(my_plan.f_hat[j]));
00092 }
00093
00094 fclose(fin);
00095 fclose(fout_real);
00096 fclose(fout_imag);
00097
00098 nfft_finalize(&my_plan);
00099 }
00100
00101
00102 int main(int argc, char **argv)
00103 {
00104 if (argc <= 5) {
00105 printf("usage: ./reconstruct_data_gridding FILENAME N M ITER WEIGHTS\n");
00106 return 1;
00107 }
00108 reconstruct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[5]));
00109
00110 return 1;
00111 }
00112