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