00001
00019 #include <stdio.h>
00020 #include <math.h>
00021 #include <stdlib.h>
00022 #include <string.h>
00023
00024 #include "util.h"
00025 #include "nfft3.h"
00026
00028
00029 #define KERNEL(r) (1.0-fabs((double)(r))/((double)R/2))
00030
00034 int polar_grid(int T, int R, double *x, double *w)
00035 {
00036 int t, r;
00037 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00038
00039 for(t=-T/2; t<T/2; t++)
00040 {
00041 for(r=-R/2; r<R/2; r++)
00042 {
00043 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
00044 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
00045 if (r==0)
00046 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00047 else
00048 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00049 }
00050 }
00051
00052 return 0;
00053 }
00054
00058 int linogram_grid(int T, int R, double *x, double *w)
00059 {
00060 int t, r;
00061 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00062
00063 for(t=-T/2; t<T/2; t++)
00064 {
00065 for(r=-R/2; r<R/2; r++)
00066 {
00067 if(t<0)
00068 {
00069 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
00070 x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
00071 }
00072 else
00073 {
00074 x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
00075 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
00076 }
00077 if (r==0)
00078 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00079 else
00080 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00081 }
00082 }
00083
00084 return 0;
00085 }
00086
00091 int Inverse_Radon_trafo(int (*gridfcn)(), int T, int R, double *Rf, int NN, double *f, int max_i)
00092 {
00093 int j,k;
00094 nfft_plan my_nfft_plan;
00095 infft_plan my_infft_plan;
00097 fftw_complex *fft;
00098 fftw_plan my_fftw_plan;
00100 int t,r;
00101 double *x, *w;
00102 int l;
00104 int N[2],n[2];
00105 int M=T*R;
00106
00107 N[0]=NN; n[0]=2*N[0];
00108 N[1]=NN; n[1]=2*N[1];
00109
00110 fft = (fftw_complex *)fftw_malloc(R*sizeof(fftw_complex));
00111 my_fftw_plan = fftw_plan_dft_1d(R,fft,fft,FFTW_FORWARD,FFTW_MEASURE);
00112
00113 x = (double *)malloc(2*T*R*(sizeof(double)));
00114 if (x==NULL)
00115 return -1;
00116
00117 w = (double *)malloc(T*R*(sizeof(double)));
00118 if (w==NULL)
00119 return -1;
00120
00122 nfft_init_guru(&my_nfft_plan, 2, N, M, n, 4,
00123 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00124 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00125
00127 infft_init_advanced(&my_infft_plan,&my_nfft_plan, CGNR | PRECOMPUTE_WEIGHT);
00128
00130 gridfcn(T,R,x,w);
00131 for(j=0;j<my_nfft_plan.M_total;j++)
00132 {
00133 my_nfft_plan.x[2*j+0] = x[2*j+0];
00134 my_nfft_plan.x[2*j+1] = x[2*j+1];
00135 if (j%R)
00136 my_infft_plan.w[j] = w[j];
00137 else
00138 my_infft_plan.w[j] = 0.0;
00139 }
00140
00142 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00143 nfft_precompute_lin_psi(&my_nfft_plan);
00144
00145 if(my_nfft_plan.nfft_flags & PRE_PSI)
00146 nfft_precompute_psi(&my_nfft_plan);
00147
00148 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00149 nfft_precompute_full_psi(&my_nfft_plan);
00150
00152 for(t=0; t<T; t++)
00153 {
00154
00155
00156
00157
00158
00159
00160 for(r=0; r<R; r++)
00161 fft[r] = Rf[t*R+r] + I*0.0;
00162
00163 nfft_fftshift_complex(fft, 1, &R);
00164 fftw_execute(my_fftw_plan);
00165 nfft_fftshift_complex(fft, 1, &R);
00166
00167 my_infft_plan.y[t*R] = 0.0;
00168 for(r=-R/2+1; r<R/2; r++)
00169 my_infft_plan.y[t*R+(r+R/2)] = fft[r+R/2]/KERNEL(r);
00170 }
00171
00173 for(k=0;k<my_nfft_plan.N_total;k++)
00174 my_infft_plan.f_hat_iter[k] = 0.0 + I*0.0;
00175
00177 infft_before_loop(&my_infft_plan);
00178
00179 if (max_i<1)
00180 {
00181 l=1;
00182 for(k=0;k<my_nfft_plan.N_total;k++)
00183 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00184 }
00185 else
00186 {
00187 for(l=1;l<=max_i;l++)
00188 {
00189 infft_loop_one_step(&my_infft_plan);
00190
00191 }
00192 }
00193
00194
00196 for(k=0;k<my_nfft_plan.N_total;k++)
00197 f[k] = creal(my_infft_plan.f_hat_iter[k]);
00198
00200 fftw_destroy_plan(my_fftw_plan);
00201 fftw_free(fft);
00202 infft_finalize(&my_infft_plan);
00203 nfft_finalize(&my_nfft_plan);
00204 free(x);
00205 free(w);
00206 return 0;
00207 }
00208
00211 int main(int argc,char **argv)
00212 {
00213 int (*gridfcn)();
00214 int T, R;
00215 FILE *fp;
00216 int N;
00217 double *Rf, *iRf;
00218 int max_i;
00220 if( argc!=6 )
00221 {
00222 printf("inverse_radon gridfcn N T R max_i\n");
00223 printf("\n");
00224 printf("gridfcn \"polar\" or \"linogram\" \n");
00225 printf("N image size NxN \n");
00226 printf("T number of slopes \n");
00227 printf("R number of offsets \n");
00228 printf("max_i number of iterations \n");
00229 exit(-1);
00230 }
00231
00232 if (strcmp(argv[1],"polar") == 0)
00233 gridfcn = polar_grid;
00234 else
00235 gridfcn = linogram_grid;
00236
00237 N = atoi(argv[2]);
00238 T = atoi(argv[3]);
00239 R = atoi(argv[4]);
00240
00241 max_i = atoi(argv[5]);
00242
00243 Rf = (double *)malloc(T*R*(sizeof(double)));
00244 iRf = (double *)malloc(N*N*(sizeof(double)));
00245
00247 fp=fopen("sinogram_data.bin","rb");
00248 if (fp==NULL)
00249 return(-1);
00250 fread(Rf,sizeof(double),T*R,fp);
00251 fclose(fp);
00252
00254 Inverse_Radon_trafo(gridfcn,T,R,Rf,N,iRf,max_i);
00255
00257 fp=fopen("output_data.bin","wb+");
00258 if (fp==NULL)
00259 return(-1);
00260 fwrite(iRf,sizeof(double),N*N,fp);
00261 fclose(fp);
00262
00264 free(Rf);
00265 free(iRf);
00266
00267 return EXIT_SUCCESS;
00268 }