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