NFFT Logo 3.2.2
radon.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 /* $Id: radon.c 3896 2012-10-10 12:19:26Z tovo $ */
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 /*#define KERNEL(r) 1.0 */
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 /*    for(r=0; r<R/2; r++)
00190       Rf[t*R+(r+R/2)] = creal(cexp(-I*PI*r)*fft[r]);
00191     for(r=0; r<R/2; r++)
00192       Rf[t*R+r] = creal(cexp(-I*PI*r)*fft[r+R/2]);
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   /*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409