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