00001 #include <stdlib.h>
00002 #include <math.h>
00003 #include "util.h"
00004 #include "nfft3.h"
00005
00015 void reconstruct(char* filename, int N, int M, int weight)
00016 {
00017 int j;
00018 double weights;
00019 double real,imag;
00020 nfft_plan my_plan;
00021 FILE* fin;
00022 FILE* fweight;
00023 FILE *fout_real;
00024 FILE *fout_imag;
00025 int my_N[2],my_n[2];
00026 int flags = PRE_PHI_HUT| PRE_PSI |MALLOC_X| MALLOC_F_HAT|
00027 MALLOC_F| FFTW_INIT| FFT_OUT_OF_PLACE|
00028 FFTW_MEASURE| FFTW_DESTROY_INPUT;
00029
00030
00031 my_N[0]=N; my_n[0]=ceil(N*1.2);
00032 my_N[1]=N; my_n[1]=ceil(N*1.2);
00033 nfft_init_guru(&my_plan, 2, my_N, M, my_n, 6,flags,
00034 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00035
00036 fin=fopen(filename,"r");
00037
00038 fweight=fopen("weights.dat","r");
00039 for(j=0;j<my_plan.M_total;j++)
00040 {
00041 fscanf(fweight,"%le ",&weights);
00042 fscanf(fin,"%le %le %le %le",&my_plan.x[2*j+0],&my_plan.x[2*j+1],&real,&imag);
00043 my_plan.f[j] = real + I*imag;
00044 if (weight)
00045 my_plan.f[j] = my_plan.f[j] * weights;
00046 }
00047 fclose(fweight);
00048
00049
00050 if(my_plan.nfft_flags & PRE_PSI)
00051 nfft_precompute_psi(&my_plan);
00052
00053
00054 if(my_plan.nfft_flags & PRE_FULL_PSI)
00055 nfft_precompute_full_psi(&my_plan);
00056
00057
00058
00059 nfft_adjoint(&my_plan);
00060
00061 fout_real=fopen("output_real.dat","w");
00062 fout_imag=fopen("output_imag.dat","w");
00063
00064 for (j=0;j<N*N;j++) {
00065 fprintf(fout_real,"%le ",creal(my_plan.f_hat[j]));
00066 fprintf(fout_imag,"%le ",cimag(my_plan.f_hat[j]));
00067 }
00068
00069 fclose(fin);
00070 fclose(fout_real);
00071 fclose(fout_imag);
00072
00073 nfft_finalize(&my_plan);
00074 }
00075
00076
00077 int main(int argc, char **argv)
00078 {
00079 if (argc <= 5) {
00080 printf("usage: ./reconstruct_data_gridding FILENAME N M ITER WEIGHTS\n");
00081 return 1;
00082 }
00083 reconstruct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[5]));
00084
00085 return 1;
00086 }
00087