inverse_radon.c

Go to the documentation of this file.

/*
 * Copyright (c) 2002, 2017 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.
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <complex.h>

#define NFFT_PRECISION_DOUBLE

#include "nfft3mp.h"

/*#define KERNEL(r) 1.0 */
#define KERNEL(r) (NFFT_K(1.0)-NFFT_M(fabs)((NFFT_R)(r))/((NFFT_R)S/2))

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

  for (t = -T / 2; t < T / 2; t++)
  {
    for (r = -S / 2; r < S / 2; r++)
    {
      x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S) * NFFT_M(cos)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
      x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S) * NFFT_M(sin)(NFFT_KPI * (NFFT_R)(t) / (NFFT_R)(T));
      if (r == 0)
        w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
      else
        w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
    }
  }

  return 0;
}

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

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

  return 0;
}

static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, NFFT_R *Rf, int NN, NFFT_R *f,
    int max_i)
{
  int j, k; 
  NFFT(plan) my_nfft_plan; 
  SOLVER(plan_complex) my_infft_plan; 
  NFFT_C *fft; 
  FFTW(plan) my_fftw_plan; 
  int t, r; 
  NFFT_R *x, *w; 
  int l; 
  int N[2], n[2];
  int M = T * S;

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

  fft = (NFFT_C *) NFFT(malloc)((size_t)(S) * sizeof(NFFT_C));
  my_fftw_plan = FFTW(plan_dft_1d)(S, fft, fft, FFTW_FORWARD, FFTW_MEASURE);

  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
  if (x == NULL)
    return EXIT_FAILURE;

  w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
  if (w == NULL)
    return EXIT_FAILURE;

  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, S, 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 % S)
      my_infft_plan.w[j] = w[j];
    else
      my_infft_plan.w[j] = NFFT_K(0.0);
  }

  if (my_nfft_plan.flags & PRE_LIN_PSI)
    NFFT(precompute_lin_psi)(&my_nfft_plan);

  if (my_nfft_plan.flags & PRE_PSI)
    NFFT(precompute_psi)(&my_nfft_plan);

  if (my_nfft_plan.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*NFFT_KPI*r)*Rf[t*R+(r+R/2)];
     for(r=0; r<R/2; r++)
     fft[r+R/2] = cexp(I*NFFT_KPI*r)*Rf[t*R+r];
     */

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

    NFFT(fftshift_complex_int)(fft, 1, &S);
    FFTW(execute)(my_fftw_plan);
    NFFT(fftshift_complex_int)(fft, 1, &S);

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

  for (k = 0; k < my_nfft_plan.N_total; k++)
    my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(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] = NFFT_M(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, S; 
  FILE *fp;
  int N; 
  NFFT_R *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(EXIT_FAILURE);
  }

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

  N = atoi(argv[2]);
  T = atoi(argv[3]);
  S = 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 = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
  iRf = (NFFT_R *) NFFT(malloc)((size_t)(N * N) * (sizeof(NFFT_R)));

  fp = fopen("sinogram_data.bin", "rb");
  if (fp == NULL)
    return EXIT_FAILURE;
  fread(Rf, sizeof(NFFT_R), (size_t)(T * S), fp);
  fclose(fp);

  inverse_radon_trafo(gridfcn, T, S, Rf, N, iRf, max_i);

  fp = fopen("output_data.bin", "wb+");
  if (fp == NULL)
    return EXIT_FAILURE;
  fwrite(iRf, sizeof(NFFT_R), (size_t)(N * N), fp);
  fclose(fp);

  NFFT(free)(Rf);
  NFFT(free)(iRf);

  return EXIT_SUCCESS;
}