00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00040 #include <stdio.h>
00041 #include <math.h>
00042 #include <stdlib.h>
00043 #include <string.h>
00044 #include <complex.h>
00045
00046 #include "nfft3util.h"
00047 #include "nfft3.h"
00048
00050
00051 #define KERNEL(r) (1.0-fabs((double)(r))/((double)R/2))
00052
00056 int polar_grid(int T, int R, double *x, double *w)
00057 {
00058 int t, r;
00059 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00060
00061 for(t=-T/2; t<T/2; t++)
00062 {
00063 for(r=-R/2; r<R/2; r++)
00064 {
00065 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
00066 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
00067 if (r==0)
00068 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00069 else
00070 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00071 }
00072 }
00073
00074 return 0;
00075 }
00076
00080 int linogram_grid(int T, int R, double *x, double *w)
00081 {
00082 int t, r;
00083 double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00084
00085 for(t=-T/2; t<T/2; t++)
00086 {
00087 for(r=-R/2; r<R/2; r++)
00088 {
00089 if(t<0)
00090 {
00091 x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
00092 x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
00093 }
00094 else
00095 {
00096 x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
00097 x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
00098 }
00099 if (r==0)
00100 w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00101 else
00102 w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00103 }
00104 }
00105
00106 return 0;
00107 }
00108
00112 int Radon_trafo(int (*gridfcn)(), int T, int R, double *f, int NN, double *Rf)
00113 {
00114 int j,k;
00115 nfft_plan my_nfft_plan;
00117 fftw_complex *fft;
00118 fftw_plan my_fftw_plan;
00120 int t,r;
00121 double *x, *w;
00123 int N[2],n[2];
00124 int M=T*R;
00125
00126 N[0]=NN; n[0]=2*N[0];
00127 N[1]=NN; n[1]=2*N[1];
00128
00129 fft = (fftw_complex *)nfft_malloc(R*sizeof(fftw_complex));
00130 my_fftw_plan = fftw_plan_dft_1d(R,fft,fft,FFTW_BACKWARD,FFTW_MEASURE);
00131
00132 x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00133 if (x==NULL)
00134 return -1;
00135
00136 w = (double *)nfft_malloc(T*R*(sizeof(double)));
00137 if (w==NULL)
00138 return -1;
00139
00141 nfft_init_guru(&my_nfft_plan, 2, N, M, n, 4,
00142 PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00143 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00144
00145
00147 gridfcn(T,R,x,w);
00148 for(j=0;j<my_nfft_plan.M_total;j++)
00149 {
00150 my_nfft_plan.x[2*j+0] = x[2*j+0];
00151 my_nfft_plan.x[2*j+1] = x[2*j+1];
00152 }
00153
00155 if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00156 nfft_precompute_lin_psi(&my_nfft_plan);
00157
00158 if(my_nfft_plan.nfft_flags & PRE_PSI)
00159 nfft_precompute_psi(&my_nfft_plan);
00160
00161 if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00162 nfft_precompute_full_psi(&my_nfft_plan);
00163
00165 for(k=0;k<my_nfft_plan.N_total;k++)
00166 my_nfft_plan.f_hat[k] = f[k] + _Complex_I*0.0;
00167
00169 nfft_trafo(&my_nfft_plan);
00170
00172 for(t=0; t<T; t++)
00173 {
00174 fft[0]=0.0;
00175 for(r=-R/2+1; r<R/2; r++)
00176 fft[r+R/2] = KERNEL(r)*my_nfft_plan.f[t*R+(r+R/2)];
00177
00178 nfft_fftshift_complex(fft, 1, &R);
00179 fftw_execute(my_fftw_plan);
00180 nfft_fftshift_complex(fft, 1, &R);
00181
00182 for(r=0; r<R; r++)
00183 Rf[t*R+r] = creal(fft[r])/R;
00184
00185
00186
00187
00188
00189
00190 }
00191
00193 fftw_destroy_plan(my_fftw_plan);
00194 nfft_free(fft);
00195 nfft_finalize(&my_nfft_plan);
00196 nfft_free(x);
00197 nfft_free(w);
00198 return 0;
00199 }
00200
00203 int main(int argc,char **argv)
00204 {
00205 int (*gridfcn)();
00206 int T, R;
00207 FILE *fp;
00208 int N;
00209 double *f, *Rf;
00210
00211 if( argc!=5 )
00212 {
00213 printf("radon gridfcn N T R\n");
00214 printf("\n");
00215 printf("gridfcn \"polar\" or \"linogram\" \n");
00216 printf("N image size NxN \n");
00217 printf("T number of slopes \n");
00218 printf("R number of offsets \n");
00219 exit(-1);
00220 }
00221
00222 if (strcmp(argv[1],"polar") == 0)
00223 gridfcn = polar_grid;
00224 else
00225 gridfcn = linogram_grid;
00226
00227 N = atoi(argv[2]);
00228 T = atoi(argv[3]);
00229 R = atoi(argv[4]);
00230
00231
00232 f = (double *)nfft_malloc(N*N*(sizeof(double)));
00233 Rf = (double *)nfft_malloc(T*R*(sizeof(double)));
00234
00236 fp=fopen("input_data.bin","rb");
00237 if (fp==NULL)
00238 return(-1);
00239 fread(f,sizeof(double),N*N,fp);
00240 fclose(fp);
00241
00243 Radon_trafo(gridfcn,T,R,f,N,Rf);
00244
00246 fp=fopen("sinogram_data.bin","wb+");
00247 if (fp==NULL)
00248 return(-1);
00249 fwrite(Rf,sizeof(double),T*R,fp);
00250 fclose(fp);
00251
00253 nfft_free(f);
00254 nfft_free(Rf);
00255
00256 return EXIT_SUCCESS;
00257 }