00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00039 #include "config.h"
00040
00041 #include <stdio.h>
00042 #include <math.h>
00043 #include <stdlib.h>
00044 #include <string.h>
00045 #ifdef HAVE_COMPLEX_H
00046 #include <complex.h>
00047 #endif
00048
00049 #include "nfft3util.h"
00050 #include "nfft3.h"
00051
00053
00054 #define KERNEL(r) (1.0-fabs((double)(r))/((double)R/2))
00055
00059 static int polar_grid(int T, int R, double *x, double *w)
00060 {
00061 int t, r;
00062 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00063
00064 for(t=-T/2; t<T/2; t++)
00065 {
00066 for(r=-R/2; r<R/2; r++)
00067 {
00068 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
00069 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
00070 if (r==0)
00071 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00072 else
00073 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00074 }
00075 }
00076
00077 return 0;
00078 }
00079
00083 static int linogram_grid(int T, int R, double *x, double *w)
00084 {
00085 int t, r;
00086 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00087
00088 for(t=-T/2; t<T/2; t++)
00089 {
00090 for(r=-R/2; r<R/2; r++)
00091 {
00092 if(t<0)
00093 {
00094 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
00095 x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
00096 }
00097 else
00098 {
00099 x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
00100 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
00101 }
00102 if (r==0)
00103 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00104 else
00105 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00106 }
00107 }
00108
00109 return 0;
00110 }
00111
00116 int Inverse_Radon_trafo(int (*gridfcn)(), int T, int R, double *Rf, int NN, double *f, int max_i)
00117 {
00118 int j,k;
00119 nfft_plan my_nfft_plan;
00120 solver_plan_complex my_infft_plan;
00122 fftw_complex *fft;
00123 fftw_plan my_fftw_plan;
00125 int t,r;
00126 double *x, *w;
00127 int l;
00129 int N[2],n[2];
00130 int M=T*R;
00131
00132 N[0]=NN; n[0]=2*N[0];
00133 N[1]=NN; n[1]=2*N[1];
00134
00135 fft = (fftw_complex *)nfft_malloc(R*sizeof(fftw_complex));
00136 my_fftw_plan = fftw_plan_dft_1d(R,fft,fft,FFTW_FORWARD,FFTW_MEASURE);
00137
00138 x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00139 if (x==NULL)
00140 return -1;
00141
00142 w = (double *)nfft_malloc(T*R*(sizeof(double)));
00143 if (w==NULL)
00144 return -1;
00145
00147 nfft_init_guru(&my_nfft_plan, 2, N, M, n, 4,
00148 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00149 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00150
00152 solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
00153
00155 gridfcn(T,R,x,w);
00156 for(j=0;j<my_nfft_plan.M_total;j++)
00157 {
00158 my_nfft_plan.x[2*j+0] = x[2*j+0];
00159 my_nfft_plan.x[2*j+1] = x[2*j+1];
00160 if (j%R)
00161 my_infft_plan.w[j] = w[j];
00162 else
00163 my_infft_plan.w[j] = 0.0;
00164 }
00165
00167 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00168 nfft_precompute_lin_psi(&my_nfft_plan);
00169
00170 if(my_nfft_plan.nfft_flags & PRE_PSI)
00171 nfft_precompute_psi(&my_nfft_plan);
00172
00173 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00174 nfft_precompute_full_psi(&my_nfft_plan);
00175
00177 for(t=0; t<T; t++)
00178 {
00179
00180
00181
00182
00183
00184
00185 for(r=0; r<R; r++)
00186 fft[r] = Rf[t*R+r] + _Complex_I*0.0;
00187
00188 nfft_fftshift_complex(fft, 1, &R);
00189 fftw_execute(my_fftw_plan);
00190 nfft_fftshift_complex(fft, 1, &R);
00191
00192 my_infft_plan.y[t*R] = 0.0;
00193 for(r=-R/2+1; r<R/2; r++)
00194 my_infft_plan.y[t*R+(r+R/2)] = fft[r+R/2]/KERNEL(r);
00195 }
00196
00198 for(k=0;k<my_nfft_plan.N_total;k++)
00199 my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
00200
00202 solver_before_loop_complex(&my_infft_plan);
00203
00204 if (max_i<1)
00205 {
00206 l=1;
00207 for(k=0;k<my_nfft_plan.N_total;k++)
00208 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00209 }
00210 else
00211 {
00212 for(l=1;l<=max_i;l++)
00213 {
00214 solver_loop_one_step_complex(&my_infft_plan);
00215
00216 }
00217 }
00218
00219
00221 for(k=0;k<my_nfft_plan.N_total;k++)
00222 f[k] = creal(my_infft_plan.f_hat_iter[k]);
00223
00225 fftw_destroy_plan(my_fftw_plan);
00226 nfft_free(fft);
00227 solver_finalize_complex(&my_infft_plan);
00228 nfft_finalize(&my_nfft_plan);
00229 nfft_free(x);
00230 nfft_free(w);
00231 return 0;
00232 }
00233
00236 int main(int argc,char **argv)
00237 {
00238 int (*gridfcn)();
00239 int T, R;
00240 FILE *fp;
00241 int N;
00242 double *Rf, *iRf;
00243 int max_i;
00245 if( argc!=6 )
00246 {
00247 printf("inverse_radon gridfcn N T R max_i\n");
00248 printf("\n");
00249 printf("gridfcn \"polar\" or \"linogram\" \n");
00250 printf("N image size NxN \n");
00251 printf("T number of slopes \n");
00252 printf("R number of offsets \n");
00253 printf("max_i number of iterations \n");
00254 exit(-1);
00255 }
00256
00257 if (strcmp(argv[1],"polar") == 0)
00258 gridfcn = polar_grid;
00259 else
00260 gridfcn = linogram_grid;
00261
00262 N = atoi(argv[2]);
00263 T = atoi(argv[3]);
00264 R = atoi(argv[4]);
00265
00266 max_i = atoi(argv[5]);
00267
00268 Rf = (double *)nfft_malloc(T*R*(sizeof(double)));
00269 iRf = (double *)nfft_malloc(N*N*(sizeof(double)));
00270
00272 fp=fopen("sinogram_data.bin","rb");
00273 if (fp==NULL)
00274 return(-1);
00275 fread(Rf,sizeof(double),T*R,fp);
00276 fclose(fp);
00277
00279 Inverse_Radon_trafo(gridfcn,T,R,Rf,N,iRf,max_i);
00280
00282 fp=fopen("output_data.bin","wb+");
00283 if (fp==NULL)
00284 return(-1);
00285 fwrite(iRf,sizeof(double),N*N,fp);
00286 fclose(fp);
00287
00289 nfft_free(Rf);
00290 nfft_free(iRf);
00291
00292 return EXIT_SUCCESS;
00293 }