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 Z, int weight ,fftw_complex *mem)
00038 {
00039 int j,k,z;
00040 double weights;
00041 double tmp;
00042 double real,imag;
00043 nfft_plan my_plan;
00044 int my_N[2],my_n[2];
00045 FILE* fin;
00046 FILE* fweight;
00047
00048
00049 my_N[0]=N; my_n[0]=ceil(N*1.2);
00050 my_N[1]=N; my_n[1]=ceil(N*1.2);
00051 nfft_init_guru(&my_plan, 2, my_N, M/Z, my_n, 6, PRE_PHI_HUT| PRE_PSI|
00052 MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00053 FFTW_INIT| FFT_OUT_OF_PLACE,
00054 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00055
00056
00057 if(my_plan.nfft_flags & PRE_LIN_PSI)
00058 nfft_precompute_lin_psi(&my_plan);
00059
00060 fin=fopen(filename,"r");
00061
00062 for(z=0;z<Z;z++) {
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 %le",
00068 &my_plan.x[2*j+0],&my_plan.x[2*j+1],&tmp,&real,&imag);
00069 my_plan.f[j] = real + _Complex_I*imag;
00070 if(weight)
00071 my_plan.f[j] = my_plan.f[j] * weights;
00072 }
00073 fclose(fweight);
00074
00075
00076 if(z==0 && my_plan.nfft_flags & PRE_PSI)
00077 nfft_precompute_psi(&my_plan);
00078
00079
00080 if(z==0 && my_plan.nfft_flags & PRE_FULL_PSI)
00081 nfft_precompute_full_psi(&my_plan);
00082
00083
00084 nfft_adjoint(&my_plan);
00085
00086 for(k=0;k<my_plan.N_total;k++) {
00087
00088
00089 mem[(Z*N*N/2+z*N*N+ k)%(Z*N*N)] = my_plan.f_hat[k];
00090 }
00091 }
00092 fclose(fin);
00093
00094 nfft_finalize(&my_plan);
00095 }
00096
00101 void print(int N,int M,int Z, fftw_complex *mem)
00102 {
00103 int i,j;
00104 FILE* fout_real;
00105 FILE* fout_imag;
00106 fout_real=fopen("output_real.dat","w");
00107 fout_imag=fopen("output_imag.dat","w");
00108
00109 for(i=0;i<Z;i++) {
00110 for (j=0;j<N*N;j++) {
00111 fprintf(fout_real,"%le ",creal(mem[(Z*N*N/2+i*N*N+ j)%(Z*N*N)]) /Z);
00112 fprintf(fout_imag,"%le ",cimag(mem[(Z*N*N/2+i*N*N+ j)%(Z*N*N)]) /Z);
00113 }
00114 fprintf(fout_real,"\n");
00115 fprintf(fout_imag,"\n");
00116 }
00117
00118 fclose(fout_real);
00119 fclose(fout_imag);
00120 }
00121
00122
00123 int main(int argc, char **argv)
00124 {
00125 fftw_complex *mem;
00126 fftw_plan plan;
00127 int N,M,Z;
00128
00129 if (argc <= 6) {
00130 printf("usage: ./reconstruct_data_gridding FILENAME N M Z ITER WEIGHTS\n");
00131 return 1;
00132 }
00133
00134 N=atoi(argv[2]);
00135 M=atoi(argv[3]);
00136 Z=atoi(argv[4]);
00137
00138
00139
00140 mem = (fftw_complex*) nfft_malloc(sizeof(fftw_complex) * atoi(argv[2]) * atoi(argv[2]) * atoi(argv[4]));
00141
00142
00143 plan = fftw_plan_many_dft(1, &Z, N*N,
00144 mem, NULL,
00145 N*N, 1,
00146 mem, NULL,
00147 N*N,1 ,
00148 FFTW_BACKWARD, FFTW_MEASURE);
00149
00150
00151 reconstruct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[4]),atoi(argv[6]),mem);
00152
00153
00154 fftw_execute(plan);
00155
00156
00157 print(N,M,Z, mem);
00158
00159
00160 nfft_free(mem);
00161
00162 return 1;
00163 }
00164