00001
00020 #include <stdio.h>
00021 #include <math.h>
00022 #include <stdlib.h>
00023 #include <string.h>
00024
00025 #include "util.h"
00026 #include "nfft3.h"
00027
00029
00030 #define KERNEL(r) (1.0-fabs((double)(r))/((double)R/2))
00031
00035 int polar_grid(int T, int R, double *x, double *w)
00036 {
00037 int t, r;
00038 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00039
00040 for(t=-T/2; t<T/2; t++)
00041 {
00042 for(r=-R/2; r<R/2; r++)
00043 {
00044 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
00045 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
00046 if (r==0)
00047 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00048 else
00049 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00050 }
00051 }
00052
00053 return 0;
00054 }
00055
00059 int linogram_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 if(t<0)
00069 {
00070 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
00071 x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
00072 }
00073 else
00074 {
00075 x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
00076 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
00077 }
00078 if (r==0)
00079 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00080 else
00081 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00082 }
00083 }
00084
00085 return 0;
00086 }
00087
00091 int Radon_trafo(int (*gridfcn)(), int T, int R, double *f, int NN, double *Rf)
00092 {
00093 int j,k;
00094 nfft_plan my_nfft_plan;
00096 fftw_complex *fft;
00097 fftw_plan my_fftw_plan;
00099 int t,r;
00100 double *x, *w;
00102 int N[2],n[2];
00103 int M=T*R;
00104
00105 N[0]=NN; n[0]=2*N[0];
00106 N[1]=NN; n[1]=2*N[1];
00107
00108 fft = (fftw_complex *)fftw_malloc(R*sizeof(fftw_complex));
00109 my_fftw_plan = fftw_plan_dft_1d(R,fft,fft,FFTW_BACKWARD,FFTW_MEASURE);
00110
00111 x = (double *)malloc(2*T*R*(sizeof(double)));
00112 if (x==NULL)
00113 return -1;
00114
00115 w = (double *)malloc(T*R*(sizeof(double)));
00116 if (w==NULL)
00117 return -1;
00118
00120 nfft_init_guru(&my_nfft_plan, 2, N, M, n, 4,
00121 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00122 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00123
00124
00126 gridfcn(T,R,x,w);
00127 for(j=0;j<my_nfft_plan.M_total;j++)
00128 {
00129 my_nfft_plan.x[2*j+0] = x[2*j+0];
00130 my_nfft_plan.x[2*j+1] = x[2*j+1];
00131 }
00132
00134 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00135 nfft_precompute_lin_psi(&my_nfft_plan);
00136
00137 if(my_nfft_plan.nfft_flags & PRE_PSI)
00138 nfft_precompute_psi(&my_nfft_plan);
00139
00140 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00141 nfft_precompute_full_psi(&my_nfft_plan);
00142
00144 for(k=0;k<my_nfft_plan.N_total;k++)
00145 my_nfft_plan.f_hat[k] = f[k] + I*0.0;
00146
00148 nfft_trafo(&my_nfft_plan);
00149
00151 for(t=0; t<T; t++)
00152 {
00153 fft[0]=0.0;
00154 for(r=-R/2+1; r<R/2; r++)
00155 fft[r+R/2] = KERNEL(r)*my_nfft_plan.f[t*R+(r+R/2)];
00156
00157 nfft_fftshift_complex(fft, 1, &R);
00158 fftw_execute(my_fftw_plan);
00159 nfft_fftshift_complex(fft, 1, &R);
00160
00161 for(r=0; r<R; r++)
00162 Rf[t*R+r] = creal(fft[r])/R;
00163
00164
00165
00166
00167
00168
00169 }
00170
00172 fftw_destroy_plan(my_fftw_plan);
00173 fftw_free(fft);
00174 nfft_finalize(&my_nfft_plan);
00175 free(x);
00176 free(w);
00177 return 0;
00178 }
00179
00182 int main(int argc,char **argv)
00183 {
00184 int (*gridfcn)();
00185 int T, R;
00186 FILE *fp;
00187 int N;
00188 double *f, *Rf;
00189
00190 if( argc!=5 )
00191 {
00192 printf("radon gridfcn N T R\n");
00193 printf("\n");
00194 printf("gridfcn \"polar\" or \"linogram\" \n");
00195 printf("N image size NxN \n");
00196 printf("T number of slopes \n");
00197 printf("R number of offsets \n");
00198 exit(-1);
00199 }
00200
00201 if (strcmp(argv[1],"polar") == 0)
00202 gridfcn = polar_grid;
00203 else
00204 gridfcn = linogram_grid;
00205
00206 N = atoi(argv[2]);
00207 T = atoi(argv[3]);
00208 R = atoi(argv[4]);
00209
00210
00211 f = (double *)malloc(N*N*(sizeof(double)));
00212 Rf = (double *)malloc(T*R*(sizeof(double)));
00213
00215 fp=fopen("input_data.bin","rb");
00216 if (fp==NULL)
00217 return(-1);
00218 fread(f,sizeof(double),N*N,fp);
00219 fclose(fp);
00220
00222 Radon_trafo(gridfcn,T,R,f,N,Rf);
00223
00225 fp=fopen("sinogram_data.bin","wb+");
00226 if (fp==NULL)
00227 return(-1);
00228 fwrite(Rf,sizeof(double),T*R,fp);
00229 fclose(fp);
00230
00232 free(f);
00233 free(Rf);
00234
00235 return EXIT_SUCCESS;
00236 }