NFFT Logo 3.2.2
inverse_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: inverse_radon.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 
00039 #include "config.h"
00040 
00041 #include <stdio.h>
00042 #include <math.h>
00043 #include <stdlib.h>
00044 #include <string.h>
00045 #ifdef HAVE_COMPLEX_H
00046 #include <complex.h>
00047 #endif
00048 
00049 #include "nfft3util.h"
00050 #include "nfft3.h"
00051 
00053 /*#define KERNEL(r) 1.0 */
00054 #define KERNEL(r) (1.0-fabs((double)(r))/((double)R/2))
00055 
00059 static int polar_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       x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
00069       x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
00070       if (r==0)
00071         w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00072       else
00073         w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00074     }
00075   }
00076 
00077   return 0;
00078 }
00079 
00083 static int linogram_grid(int T, int R, double *x, double *w)
00084 {
00085   int t, r;
00086   double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);
00087 
00088   for(t=-T/2; t<T/2; t++)
00089   {
00090     for(r=-R/2; r<R/2; r++)
00091     {
00092       if(t<0)
00093       {
00094         x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
00095         x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
00096       }
00097       else
00098       {
00099         x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
00100         x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
00101       }
00102       if (r==0)
00103         w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
00104       else
00105         w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
00106     }
00107   }
00108 
00109   return 0;
00110 }
00111 
00116 int Inverse_Radon_trafo(int (*gridfcn)(), int T, int R, double *Rf, int NN, double *f, int max_i)
00117 {
00118   int j,k;                              
00119   nfft_plan my_nfft_plan;               
00120   solver_plan_complex my_infft_plan;             
00122   fftw_complex *fft;                    
00123   fftw_plan my_fftw_plan;               
00125   int t,r;                              
00126   double *x, *w;                        
00127   int l;                                
00129   int N[2],n[2];
00130   int M=T*R;
00131 
00132   N[0]=NN; n[0]=2*N[0];
00133   N[1]=NN; n[1]=2*N[1];
00134 
00135   fft = (fftw_complex *)nfft_malloc(R*sizeof(fftw_complex));
00136   my_fftw_plan = fftw_plan_dft_1d(R,fft,fft,FFTW_FORWARD,FFTW_MEASURE);
00137 
00138   x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
00139   if (x==NULL)
00140     return -1;
00141 
00142   w = (double *)nfft_malloc(T*R*(sizeof(double)));
00143   if (w==NULL)
00144     return -1;
00145 
00147   nfft_init_guru(&my_nfft_plan, 2, N, M, n, 4,
00148                   PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
00149                   FFTW_MEASURE| FFTW_DESTROY_INPUT);
00150 
00152   solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
00153 
00155   gridfcn(T,R,x,w);
00156   for(j=0;j<my_nfft_plan.M_total;j++)
00157   {
00158     my_nfft_plan.x[2*j+0] = x[2*j+0];
00159     my_nfft_plan.x[2*j+1] = x[2*j+1];
00160     if (j%R)
00161       my_infft_plan.w[j]    = w[j];
00162     else
00163       my_infft_plan.w[j]    = 0.0;
00164   }
00165 
00167   if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
00168     nfft_precompute_lin_psi(&my_nfft_plan);
00169 
00170   if(my_nfft_plan.nfft_flags & PRE_PSI)
00171     nfft_precompute_psi(&my_nfft_plan);
00172 
00173   if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
00174     nfft_precompute_full_psi(&my_nfft_plan);
00175 
00177   for(t=0; t<T; t++)
00178   {
00179 /*    for(r=0; r<R/2; r++)
00180        fft[r] = cexp(I*PI*r)*Rf[t*R+(r+R/2)];
00181       for(r=0; r<R/2; r++)
00182        fft[r+R/2] = cexp(I*PI*r)*Rf[t*R+r];
00183  */
00184 
00185     for(r=0; r<R; r++)
00186       fft[r] = Rf[t*R+r] + _Complex_I*0.0;
00187 
00188     nfft_fftshift_complex(fft, 1, &R);
00189     fftw_execute(my_fftw_plan);
00190     nfft_fftshift_complex(fft, 1, &R);
00191 
00192     my_infft_plan.y[t*R] = 0.0;
00193     for(r=-R/2+1; r<R/2; r++)
00194       my_infft_plan.y[t*R+(r+R/2)] = fft[r+R/2]/KERNEL(r);
00195   }
00196 
00198   for(k=0;k<my_nfft_plan.N_total;k++)
00199     my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;
00200 
00202   solver_before_loop_complex(&my_infft_plan);
00203 
00204   if (max_i<1)
00205   {
00206     l=1;
00207     for(k=0;k<my_nfft_plan.N_total;k++)
00208       my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00209   }
00210   else
00211   {
00212     for(l=1;l<=max_i;l++)
00213     {
00214       solver_loop_one_step_complex(&my_infft_plan);
00215       /*if (sqrt(my_infft_plan.dot_r_iter)<=1e-12) break;*/
00216     }
00217   }
00218   /*printf("after %d iteration(s): weighted 2-norm of original residual vector = %g\n",l-1,sqrt(my_infft_plan.dot_r_iter));*/
00219 
00221   for(k=0;k<my_nfft_plan.N_total;k++)
00222     f[k] = creal(my_infft_plan.f_hat_iter[k]);
00223 
00225   fftw_destroy_plan(my_fftw_plan);
00226   nfft_free(fft);
00227   solver_finalize_complex(&my_infft_plan);
00228   nfft_finalize(&my_nfft_plan);
00229   nfft_free(x);
00230   nfft_free(w);
00231   return 0;
00232 }
00233 
00236 int main(int argc,char **argv)
00237 {
00238   int (*gridfcn)();                     
00239   int T, R;                             
00240   FILE *fp;
00241   int N;                                
00242   double *Rf, *iRf;
00243   int max_i;                            
00245   if( argc!=6 )
00246   {
00247     printf("inverse_radon gridfcn N T R max_i\n");
00248     printf("\n");
00249     printf("gridfcn    \"polar\" or \"linogram\" \n");
00250     printf("N          image size NxN            \n");
00251     printf("T          number of slopes          \n");
00252     printf("R          number of offsets         \n");
00253     printf("max_i      number of iterations      \n");
00254     exit(-1);
00255   }
00256 
00257   if (strcmp(argv[1],"polar") == 0)
00258     gridfcn = polar_grid;
00259   else
00260     gridfcn = linogram_grid;
00261 
00262   N = atoi(argv[2]);
00263   T = atoi(argv[3]);
00264   R = atoi(argv[4]);
00265   /*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
00266   max_i = atoi(argv[5]);
00267 
00268   Rf  = (double *)nfft_malloc(T*R*(sizeof(double)));
00269   iRf = (double *)nfft_malloc(N*N*(sizeof(double)));
00270 
00272   fp=fopen("sinogram_data.bin","rb");
00273   if (fp==NULL)
00274     return(-1);
00275   fread(Rf,sizeof(double),T*R,fp);
00276   fclose(fp);
00277 
00279   Inverse_Radon_trafo(gridfcn,T,R,Rf,N,iRf,max_i);
00280 
00282   fp=fopen("output_data.bin","wb+");
00283   if (fp==NULL)
00284     return(-1);
00285   fwrite(iRf,sizeof(double),N*N,fp);
00286   fclose(fp);
00287 
00289   nfft_free(Rf);
00290   nfft_free(iRf);
00291 
00292   return EXIT_SUCCESS;
00293 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409