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 <math.h>
00023 #include <stdlib.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 iteration, int weight)
00041 {
00042 int j,k,z,l;
00043 double real,imag;
00044 nfft_plan my_plan;
00045 solver_plan_complex my_iplan;
00046 FILE* fin;
00047 FILE* fout_real;
00048 FILE* fout_imag;
00049 int my_N[3],my_n[3];
00050 double epsilon=0.0000003;
00051
00052
00053 unsigned infft_flags = CGNR | PRECOMPUTE_DAMP;
00054
00055
00056
00057 my_N[0]=Z; my_n[0]=ceil(Z*1.2);
00058 my_N[1]=N; my_n[1]=ceil(N*1.2);
00059 my_N[2]=N; my_n[2]=ceil(N*1.2);
00060 nfft_init_guru(&my_plan, 3, my_N, M, my_n, 6,
00061 PRE_PHI_HUT| PRE_PSI |MALLOC_X| MALLOC_F_HAT|
00062 MALLOC_F| FFTW_INIT| FFT_OUT_OF_PLACE,
00063 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00064
00065
00066 if(my_plan.nfft_flags & PRE_LIN_PSI)
00067 nfft_precompute_lin_psi(&my_plan);
00068
00069 if (weight)
00070 infft_flags = infft_flags | PRECOMPUTE_WEIGHT;
00071
00072
00073 solver_init_advanced_complex(&my_iplan,(nfft_mv_plan_complex*)(&my_plan), infft_flags );
00074
00075
00076 if(my_iplan.flags & PRECOMPUTE_WEIGHT)
00077 {
00078 fin=fopen("weights.dat","r");
00079 for(j=0;j<M;j++)
00080 {
00081 fscanf(fin,"%le ",&my_iplan.w[j]);
00082 }
00083 fclose(fin);
00084 }
00085
00086
00087 if(my_iplan.flags & PRECOMPUTE_DAMP)
00088 {
00089 for(j=0;j<N;j++){
00090 for(k=0;k<N;k++) {
00091 for(z=0;z<N;z++) {
00092 int j2= j-N/2;
00093 int k2= k-N/2;
00094 int z2= z-N/2;
00095 double r=sqrt(j2*j2+k2*k2+z2*z2);
00096 if(r>(double) N/2)
00097 my_iplan.w_hat[z*N*N+j*N+k]=0.0;
00098 else
00099 my_iplan.w_hat[z*N*N+j*N+k]=1.0;
00100 }
00101 }
00102 }
00103 }
00104
00105
00106 fin=fopen(filename,"r");
00107
00108
00109 fout_real=fopen("output_real.dat","w");
00110 fout_imag=fopen("output_imag.dat","w");
00111
00112
00113 for(j=0;j<M;j++)
00114 {
00115 fscanf(fin,"%le %le %le %le %le ",&my_plan.x[3*j+1],&my_plan.x[3*j+2], &my_plan.x[3*j+0],
00116 &real,&imag);
00117 my_iplan.y[j] = real + _Complex_I*imag;
00118 }
00119
00120
00121 if(my_plan.nfft_flags & PRE_PSI)
00122 nfft_precompute_psi(&my_plan);
00123
00124
00125 if(my_plan.nfft_flags & PRE_FULL_PSI)
00126 nfft_precompute_full_psi(&my_plan);
00127
00128
00129 for(k=0;k<my_plan.N_total;k++)
00130 my_iplan.f_hat_iter[k]=0.0;
00131
00132
00133 solver_before_loop_complex(&my_iplan);
00134 for(l=0;l<iteration;l++)
00135 {
00136
00137 if(my_iplan.dot_r_iter<epsilon)
00138 break;
00139 fprintf(stderr,"%e, %i of %i\n",sqrt(my_iplan.dot_r_iter),
00140 l+1,iteration);
00141 solver_loop_one_step_complex(&my_iplan);
00142 }
00143
00144 for(l=0;l<Z;l++)
00145 {
00146 for(k=0;k<N*N;k++)
00147 {
00148
00149 fprintf(fout_real,"%le ",creal(my_iplan.f_hat_iter[ k+N*N*l ]));
00150 fprintf(fout_imag,"%le ",cimag(my_iplan.f_hat_iter[ k+N*N*l ]));
00151 }
00152 fprintf(fout_real,"\n");
00153 fprintf(fout_imag,"\n");
00154 }
00155
00156 fclose(fout_real);
00157 fclose(fout_imag);
00158
00159 solver_finalize_complex(&my_iplan);
00160 nfft_finalize(&my_plan);
00161 }
00162
00163 int main(int argc, char **argv)
00164 {
00165 if (argc <= 6) {
00166 printf("usage: ./reconstruct3D FILENAME N M Z ITER WEIGHTS\n");
00167 return 1;
00168 }
00169
00170 reconstruct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[4]),atoi(argv[5]),atoi(argv[6]));
00171 return 1;
00172 }
00173