NFFT  3.3.1
reconstruct_data_inh_nnfft.c
00001 /*
00002  * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 #include <stdlib.h>
00019 #include <math.h>
00020 #include <limits.h>
00021 #include <complex.h>
00022 
00023 #include "nfft3.h"
00024 
00025 #ifndef MAX
00026 #define MAX(a,b) (((a)>(b))?(a):(b))
00027 #endif
00028 
00038 static void reconstruct(char* filename,int N,int M,int iteration, int weight)
00039 {
00040   int j,k,l;                    /* some variables  */
00041   nnfft_plan my_plan;            /* plan for the two dimensional nfft  */
00042   solver_plan_complex my_iplan;          /* plan for the two dimensional infft */
00043   FILE* fin;                    /* input file                         */
00044   FILE* finh;
00045   FILE* ftime;
00046   FILE* fout_real;              /* output file                        */
00047   FILE* fout_imag;              /* output file                        */
00048   int my_N[3],my_n[3];          /* to init the nfft */
00049   double t0, t1;
00050   double t,epsilon=0.0000003;     /* epsilon is a the break criterium for
00051                                    the iteration */
00052   unsigned infft_flags = CGNR | PRECOMPUTE_DAMP; /* flags for the infft*/
00053   double time,min_time,max_time,min_inh,max_inh;
00054   double real,imag;
00055   double *w;
00056 
00057   double Ts;
00058   double W;
00059   int N3;
00060   int m=2;
00061   double sigma = 1.25;
00062 
00063   w = (double*)nfft_malloc(N*N*sizeof(double));
00064 
00065   ftime=fopen("readout_time.dat","r");
00066   finh=fopen("inh.dat","r");
00067 
00068   min_time=INT_MAX; max_time=INT_MIN;
00069   for(j=0;j<M;j++)
00070   {
00071     fscanf(ftime,"%le ",&time);
00072     if(time<min_time)
00073       min_time = time;
00074     if(time>max_time)
00075       max_time = time;
00076   }
00077 
00078   fclose(ftime);
00079 
00080   Ts=(min_time+max_time)/2.0;
00081 
00082   min_inh=INT_MAX; max_inh=INT_MIN;
00083   for(j=0;j<N*N;j++)
00084   {
00085     fscanf(finh,"%le ",&w[j]);
00086     if(w[j]<min_inh)
00087       min_inh = w[j];
00088     if(w[j]>max_inh)
00089       max_inh = w[j];
00090   }
00091   fclose(finh);
00092 
00093   N3=ceil((MAX(fabs(min_inh),fabs(max_inh))*(max_time-min_time)/2.0)*4);
00094 
00095 
00096   W=MAX(fabs(min_inh),fabs(max_inh))*2.0;
00097 
00098   fprintf(stderr,"3:  %i %e %e %e %e %e %e\n",N3,W,min_inh,max_inh,min_time,max_time,Ts);
00099 
00100   /* initialise my_plan */
00101   my_N[0]=N;my_n[0]=ceil(N*sigma);
00102   my_N[1]=N; my_n[1]=ceil(N*sigma);
00103   my_N[2]=N3; my_n[2]=ceil(N3*sigma);
00104   nnfft_init_guru(&my_plan, 3, N*N, M, my_N,my_n,m,
00105         PRE_PSI| PRE_PHI_HUT| MALLOC_X| MALLOC_V| MALLOC_F_HAT| MALLOC_F );
00106 
00107   /* precompute lin psi if set */
00108   if(my_plan.nnfft_flags & PRE_LIN_PSI)
00109     nnfft_precompute_lin_psi(&my_plan);
00110 
00111   /* set the flags for the infft*/
00112   if (weight)
00113     infft_flags = infft_flags | PRECOMPUTE_WEIGHT;
00114 
00115   /* initialise my_iplan, advanced */
00116   solver_init_advanced_complex(&my_iplan,(nfft_mv_plan_complex*)(&my_plan), infft_flags );
00117 
00118   /* get the weights */
00119   if(my_iplan.flags & PRECOMPUTE_WEIGHT)
00120   {
00121     fin=fopen("weights.dat","r");
00122     for(j=0;j<my_plan.M_total;j++)
00123     {
00124         fscanf(fin,"%le ",&my_iplan.w[j]);
00125     }
00126     fclose(fin);
00127   }
00128 
00129   /* get the damping factors */
00130   if(my_iplan.flags & PRECOMPUTE_DAMP)
00131   {
00132     for(j=0;j<N;j++){
00133       for(k=0;k<N;k++) {
00134         int j2= j-N/2;
00135         int k2= k-N/2;
00136         double r=sqrt(j2*j2+k2*k2);
00137         if(r>(double) N/2)
00138           my_iplan.w_hat[j*N+k]=0.0;
00139         else
00140           my_iplan.w_hat[j*N+k]=1.0;
00141       }
00142     }
00143   }
00144 
00145   /* open the input file */
00146   fin=fopen(filename,"r");
00147   ftime=fopen("readout_time.dat","r");
00148 
00149   for(j=0;j<my_plan.M_total;j++)
00150   {
00151     fscanf(fin,"%le %le %le %le ",&my_plan.x[3*j+0],&my_plan.x[3*j+1],&real,&imag);
00152     my_iplan.y[j]=real+ _Complex_I*imag;
00153     fscanf(ftime,"%le ",&my_plan.x[3*j+2]);
00154 
00155     my_plan.x[3*j+2] = (my_plan.x[3*j+2]-Ts)*W/N3;
00156   }
00157 
00158   for(j=0;j<N;j++)
00159     {
00160     for(l=0;l<N;l++)
00161       {
00162         my_plan.v[3*(N*j+l)+0]=(((double) j) -(((double) N)/2.0))/((double) N);
00163         my_plan.v[3*(N*j+l)+1]=(((double) l) -(((double) N)/2.0))/((double) N);
00164         my_plan.v[3*(N*j+l)+2] = w[N*j+l]/W ;
00165       }
00166     }
00167 
00168   /* precompute psi */
00169   if(my_plan.nnfft_flags & PRE_PSI) {
00170     nnfft_precompute_psi(&my_plan);
00171     if(my_plan.nnfft_flags & PRE_FULL_PSI)
00172       nnfft_precompute_full_psi(&my_plan);
00173   }
00174 
00175   if(my_plan.nnfft_flags & PRE_PHI_HUT)
00176     nnfft_precompute_phi_hut(&my_plan);
00177 
00178   /* init some guess */
00179   for(k=0;k<my_plan.N_total;k++)
00180   {
00181     my_iplan.f_hat_iter[k]=0.0;
00182   }
00183 
00184   t0 = nfft_clock_gettime_seconds();
00185 
00186   /* inverse trafo */
00187   solver_before_loop_complex(&my_iplan);
00188   for(l=0;l<iteration;l++)
00189   {
00190     /* break if dot_r_iter is smaller than epsilon*/
00191     if(my_iplan.dot_r_iter<epsilon)
00192     break;
00193     fprintf(stderr,"%e,  %i of %i\n",sqrt(my_iplan.dot_r_iter),
00194     l+1,iteration);
00195     solver_loop_one_step_complex(&my_iplan);
00196   }
00197 
00198   t1 = nfft_clock_gettime_seconds();
00199   t = t1-t0;
00200 
00201   fout_real=fopen("output_real.dat","w");
00202   fout_imag=fopen("output_imag.dat","w");
00203 
00204   for(k=0;k<my_plan.N_total;k++) {
00205 
00206     my_iplan.f_hat_iter[k]*=cexp(2.0*_Complex_I*M_PI*Ts*w[k]);
00207 
00208     fprintf(fout_real,"%le ", creal(my_iplan.f_hat_iter[k]));
00209     fprintf(fout_imag,"%le ", cimag(my_iplan.f_hat_iter[k]));
00210   }
00211 
00212 
00213   fclose(fout_real);
00214   fclose(fout_imag);
00215 
00216 
00217   /* finalize the infft */
00218   solver_finalize_complex(&my_iplan);
00219 
00220   /* finalize the nfft */
00221   nnfft_finalize(&my_plan);
00222 
00223   nfft_free(w);
00224 }
00225 
00226 int main(int argc, char **argv)
00227 {
00228   if (argc <= 5) {
00229     printf("usage: ./reconstruct_data_inh_nnfft FILENAME N M ITER WEIGHTS\n");
00230     return 1;
00231   }
00232 
00233   reconstruct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[4]),atoi(argv[5]));
00234 
00235   return 1;
00236 }
00237 /* \} */