inverse_radon.c

Go to the documentation of this file.

/*
 * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
 *
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software
 * Foundation; either version 2 of the License, or (at your option) any later
 * version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 * details.
 *
 * You should have received a copy of the GNU General Public License along with
 * this program; if not, write to the Free Software Foundation, Inc., 51
 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* $Id: inverse_radon.c 3896 2012-10-10 12:19:26Z tovo $ */

#include "config.h"

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#ifdef HAVE_COMPLEX_H
#include <complex.h>
#endif

#include "nfft3util.h"
#include "nfft3.h"

/*#define KERNEL(r) 1.0 */
#define KERNEL(r) (1.0-fabs((double)(r))/((double)R/2))

static int polar_grid(int T, int R, double *x, double *w)
{
  int t, r;
  double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);

  for(t=-T/2; t<T/2; t++)
  {
    for(r=-R/2; r<R/2; r++)
    {
      x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R*cos(PI*t/T);
      x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R*sin(PI*t/T);
      if (r==0)
        w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
      else
        w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
    }
  }

  return 0;
}

static int linogram_grid(int T, int R, double *x, double *w)
{
  int t, r;
  double W=(double)T*(((double)R/2.0)*((double)R/2.0)+1.0/4.0);

  for(t=-T/2; t<T/2; t++)
  {
    for(r=-R/2; r<R/2; r++)
    {
      if(t<0)
      {
        x[2*((t+T/2)*R+(r+R/2))+0] = (double)r/R;
        x[2*((t+T/2)*R+(r+R/2))+1] = (double)4*(t+T/4)/T*r/R;
      }
      else
      {
        x[2*((t+T/2)*R+(r+R/2))+0] = -(double)4*(t-T/4)/T*r/R;
        x[2*((t+T/2)*R+(r+R/2))+1] = (double)r/R;
      }
      if (r==0)
        w[(t+T/2)*R+(r+R/2)] = 1.0/4.0/W;
      else
        w[(t+T/2)*R+(r+R/2)] = fabs((double)r)/W;
    }
  }

  return 0;
}

int Inverse_Radon_trafo(int (*gridfcn)(), int T, int R, double *Rf, int NN, double *f, int max_i)
{
  int j,k;                              
  nfft_plan my_nfft_plan;               
  solver_plan_complex my_infft_plan;             
  fftw_complex *fft;                    
  fftw_plan my_fftw_plan;               
  int t,r;                              
  double *x, *w;                        
  int l;                                
  int N[2],n[2];
  int M=T*R;

  N[0]=NN; n[0]=2*N[0];
  N[1]=NN; n[1]=2*N[1];

  fft = (fftw_complex *)nfft_malloc(R*sizeof(fftw_complex));
  my_fftw_plan = fftw_plan_dft_1d(R,fft,fft,FFTW_FORWARD,FFTW_MEASURE);

  x = (double *)nfft_malloc(2*T*R*(sizeof(double)));
  if (x==NULL)
    return -1;

  w = (double *)nfft_malloc(T*R*(sizeof(double)));
  if (w==NULL)
    return -1;

  nfft_init_guru(&my_nfft_plan, 2, N, M, n, 4,
                  PRE_PHI_HUT| PRE_PSI| MALLOC_X | MALLOC_F_HAT| MALLOC_F| FFTW_INIT | FFT_OUT_OF_PLACE,
                  FFTW_MEASURE| FFTW_DESTROY_INPUT);

  solver_init_advanced_complex(&my_infft_plan,(nfft_mv_plan_complex*)(&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);

  gridfcn(T,R,x,w);
  for(j=0;j<my_nfft_plan.M_total;j++)
  {
    my_nfft_plan.x[2*j+0] = x[2*j+0];
    my_nfft_plan.x[2*j+1] = x[2*j+1];
    if (j%R)
      my_infft_plan.w[j]    = w[j];
    else
      my_infft_plan.w[j]    = 0.0;
  }

  if(my_nfft_plan.nfft_flags & PRE_LIN_PSI)
    nfft_precompute_lin_psi(&my_nfft_plan);

  if(my_nfft_plan.nfft_flags & PRE_PSI)
    nfft_precompute_psi(&my_nfft_plan);

  if(my_nfft_plan.nfft_flags & PRE_FULL_PSI)
    nfft_precompute_full_psi(&my_nfft_plan);

  for(t=0; t<T; t++)
  {
/*    for(r=0; r<R/2; r++)
       fft[r] = cexp(I*PI*r)*Rf[t*R+(r+R/2)];
      for(r=0; r<R/2; r++)
       fft[r+R/2] = cexp(I*PI*r)*Rf[t*R+r];
 */

    for(r=0; r<R; r++)
      fft[r] = Rf[t*R+r] + _Complex_I*0.0;

    nfft_fftshift_complex(fft, 1, &R);
    fftw_execute(my_fftw_plan);
    nfft_fftshift_complex(fft, 1, &R);

    my_infft_plan.y[t*R] = 0.0;
    for(r=-R/2+1; r<R/2; r++)
      my_infft_plan.y[t*R+(r+R/2)] = fft[r+R/2]/KERNEL(r);
  }

  for(k=0;k<my_nfft_plan.N_total;k++)
    my_infft_plan.f_hat_iter[k] = 0.0 + _Complex_I*0.0;

  solver_before_loop_complex(&my_infft_plan);

  if (max_i<1)
  {
    l=1;
    for(k=0;k<my_nfft_plan.N_total;k++)
      my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
  }
  else
  {
    for(l=1;l<=max_i;l++)
    {
      solver_loop_one_step_complex(&my_infft_plan);
      /*if (sqrt(my_infft_plan.dot_r_iter)<=1e-12) break;*/
    }
  }
  /*printf("after %d iteration(s): weighted 2-norm of original residual vector = %g\n",l-1,sqrt(my_infft_plan.dot_r_iter));*/

  for(k=0;k<my_nfft_plan.N_total;k++)
    f[k] = creal(my_infft_plan.f_hat_iter[k]);

  fftw_destroy_plan(my_fftw_plan);
  nfft_free(fft);
  solver_finalize_complex(&my_infft_plan);
  nfft_finalize(&my_nfft_plan);
  nfft_free(x);
  nfft_free(w);
  return 0;
}

int main(int argc,char **argv)
{
  int (*gridfcn)();                     
  int T, R;                             
  FILE *fp;
  int N;                                
  double *Rf, *iRf;
  int max_i;                            
  if( argc!=6 )
  {
    printf("inverse_radon gridfcn N T R max_i\n");
    printf("\n");
    printf("gridfcn    \"polar\" or \"linogram\" \n");
    printf("N          image size NxN            \n");
    printf("T          number of slopes          \n");
    printf("R          number of offsets         \n");
    printf("max_i      number of iterations      \n");
    exit(-1);
  }

  if (strcmp(argv[1],"polar") == 0)
    gridfcn = polar_grid;
  else
    gridfcn = linogram_grid;

  N = atoi(argv[2]);
  T = atoi(argv[3]);
  R = atoi(argv[4]);
  /*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
  max_i = atoi(argv[5]);

  Rf  = (double *)nfft_malloc(T*R*(sizeof(double)));
  iRf = (double *)nfft_malloc(N*N*(sizeof(double)));

  fp=fopen("sinogram_data.bin","rb");
  if (fp==NULL)
    return(-1);
  fread(Rf,sizeof(double),T*R,fp);
  fclose(fp);

  Inverse_Radon_trafo(gridfcn,T,R,Rf,N,iRf,max_i);

  fp=fopen("output_data.bin","wb+");
  if (fp==NULL)
    return(-1);
  fwrite(iRf,sizeof(double),N*N,fp);
  fclose(fp);

  nfft_free(Rf);
  nfft_free(iRf);

  return EXIT_SUCCESS;
}