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