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 Z, int weight ,fftw_complex *mem)
00016 {
00017 int j,k,z;
00018 double weights;
00019 double tmp;
00020 double real,imag;
00021 nfft_plan my_plan;
00022 int my_N[2],my_n[2];
00023 FILE* fin;
00024 FILE* fweight;
00025
00026
00027 my_N[0]=N; my_n[0]=ceil(N*1.2);
00028 my_N[1]=N; my_n[1]=ceil(N*1.2);
00029 nfft_init_guru(&my_plan, 2, my_N, M/Z, my_n, 6, PRE_PHI_HUT| PRE_PSI|
00030 MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00031 FFTW_INIT| FFT_OUT_OF_PLACE,
00032 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00033
00034
00035 if(my_plan.nfft_flags & PRE_LIN_PSI)
00036 nfft_precompute_lin_psi(&my_plan);
00037
00038 fin=fopen(filename,"r");
00039
00040 for(z=0;z<Z;z++) {
00041 fweight=fopen("weights.dat","r");
00042 for(j=0;j<my_plan.M_total;j++)
00043 {
00044 fscanf(fweight,"%le ",&weights);
00045 fscanf(fin,"%le %le %le %le %le",
00046 &my_plan.x[2*j+0],&my_plan.x[2*j+1],&tmp,&real,&imag);
00047 my_plan.f[j] = real + I*imag;
00048 if(weight)
00049 my_plan.f[j] = my_plan.f[j] * weights;
00050 }
00051 fclose(fweight);
00052
00053
00054 if(z==0 && my_plan.nfft_flags & PRE_PSI)
00055 nfft_precompute_psi(&my_plan);
00056
00057
00058 if(z==0 && my_plan.nfft_flags & PRE_FULL_PSI)
00059 nfft_precompute_full_psi(&my_plan);
00060
00061
00062 nfft_adjoint(&my_plan);
00063
00064 for(k=0;k<my_plan.N_total;k++) {
00065
00066
00067 mem[(Z*N*N/2+z*N*N+ k)%(Z*N*N)] = my_plan.f_hat[k];
00068 }
00069 }
00070 fclose(fin);
00071
00072 nfft_finalize(&my_plan);
00073 }
00074
00079 void print(int N,int M,int Z, fftw_complex *mem)
00080 {
00081 int i,j;
00082 FILE* fout_real;
00083 FILE* fout_imag;
00084 fout_real=fopen("output_real.dat","w");
00085 fout_imag=fopen("output_imag.dat","w");
00086
00087 for(i=0;i<Z;i++) {
00088 for (j=0;j<N*N;j++) {
00089 fprintf(fout_real,"%le ",creal(mem[(Z*N*N/2+i*N*N+ j)%(Z*N*N)]) /Z);
00090 fprintf(fout_imag,"%le ",cimag(mem[(Z*N*N/2+i*N*N+ j)%(Z*N*N)]) /Z);
00091 }
00092 fprintf(fout_real,"\n");
00093 fprintf(fout_imag,"\n");
00094 }
00095
00096 fclose(fout_real);
00097 fclose(fout_imag);
00098 }
00099
00100
00101 int main(int argc, char **argv)
00102 {
00103 fftw_complex *mem;
00104 fftw_plan plan;
00105 int N,M,Z;
00106
00107 if (argc <= 6) {
00108 printf("usage: ./reconstruct_data_gridding FILENAME N M Z ITER WEIGHTS\n");
00109 return 1;
00110 }
00111
00112 N=atoi(argv[2]);
00113 M=atoi(argv[3]);
00114 Z=atoi(argv[4]);
00115
00116
00117
00118 mem = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * atoi(argv[2]) * atoi(argv[2]) * atoi(argv[4]));
00119
00120
00121 plan = fftw_plan_many_dft(1, &Z, N*N,
00122 mem, NULL,
00123 N*N, 1,
00124 mem, NULL,
00125 N*N,1 ,
00126 FFTW_BACKWARD, FFTW_MEASURE);
00127
00128
00129 reconstruct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[4]),atoi(argv[6]),mem);
00130
00131
00132 fftw_execute(plan);
00133
00134
00135 print(N,M,Z, mem);
00136
00137
00138 fftw_free(mem);
00139
00140 return 1;
00141 }
00142