NFFT  3.3.1
inverse_radon.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2002, 2016 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 
00038 #include <stdio.h>
00039 #include <math.h>
00040 #include <stdlib.h>
00041 #include <string.h>
00042 #include <complex.h>
00043 
00044 #define NFFT_PRECISION_DOUBLE
00045 
00046 #include "nfft3mp.h"
00047 
00049 /*#define KERNEL(r) 1.0 */
00050 #define KERNEL(r) (NFFT_K(1.0)-NFFT_M(fabs)((NFFT_R)(r))/((NFFT_R)S/2))
00051 
00055 static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
00056 {
00057   int t, r;
00058   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));
00059 
00060   for (t = -T / 2; t < T / 2; t++)
00061   {
00062     for (r = -S / 2; r < S / 2; r++)
00063     {
00064       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));
00065       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));
00066       if (r == 0)
00067         w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
00068       else
00069         w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
00070     }
00071   }
00072 
00073   return 0;
00074 }
00075 
00079 static int linogram_grid(int T, int S, NFFT_R *x, NFFT_R *w)
00080 {
00081   int t, r;
00082   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));
00083 
00084   for (t = -T / 2; t < T / 2; t++)
00085   {
00086     for (r = -S / 2; r < S / 2; r++)
00087     {
00088       if (t < 0)
00089       {
00090         x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S);
00091         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)
00092             / (NFFT_R)(S);
00093       }
00094       else
00095       {
00096         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)
00097             * (NFFT_R)(r) / (NFFT_R)(S);
00098         x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S);
00099       }
00100       if (r == 0)
00101         w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
00102       else
00103         w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
00104     }
00105   }
00106 
00107   return 0;
00108 }
00109 
00114 static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, NFFT_R *Rf, int NN, NFFT_R *f,
00115     int max_i)
00116 {
00117   int j, k; 
00118   NFFT(plan) my_nfft_plan; 
00119   SOLVER(plan_complex) my_infft_plan; 
00121   NFFT_C *fft; 
00122   FFTW(plan) my_fftw_plan; 
00124   int t, r; 
00125   NFFT_R *x, *w; 
00126   int l; 
00128   int N[2], n[2];
00129   int M = T * S;
00130 
00131   N[0] = NN;
00132   n[0] = 2 * N[0];
00133   N[1] = NN;
00134   n[1] = 2 * N[1];
00135 
00136   fft = (NFFT_C *) NFFT(malloc)((size_t)(S) * sizeof(NFFT_C));
00137   my_fftw_plan = FFTW(plan_dft_1d)(S, fft, fft, FFTW_FORWARD, FFTW_MEASURE);
00138 
00139   x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
00140   if (x == NULL)
00141     return EXIT_FAILURE;
00142 
00143   w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
00144   if (w == NULL)
00145     return EXIT_FAILURE;
00146 
00148   NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, 4,
00149       PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
00150           | FFT_OUT_OF_PLACE,
00151       FFTW_MEASURE | FFTW_DESTROY_INPUT);
00152 
00154   SOLVER(init_advanced_complex)(&my_infft_plan,
00155       (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
00156 
00158   gridfcn(T, S, x, w);
00159   for (j = 0; j < my_nfft_plan.M_total; j++)
00160   {
00161     my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
00162     my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
00163     if (j % S)
00164       my_infft_plan.w[j] = w[j];
00165     else
00166       my_infft_plan.w[j] = NFFT_K(0.0);
00167   }
00168 
00170   if (my_nfft_plan.flags & PRE_LIN_PSI)
00171     NFFT(precompute_lin_psi)(&my_nfft_plan);
00172 
00173   if (my_nfft_plan.flags & PRE_PSI)
00174     NFFT(precompute_psi)(&my_nfft_plan);
00175 
00176   if (my_nfft_plan.flags & PRE_FULL_PSI)
00177     NFFT(precompute_full_psi)(&my_nfft_plan);
00178 
00180   for (t = 0; t < T; t++)
00181   {
00182     /*    for(r=0; r<R/2; r++)
00183      fft[r] = cexp(I*NFFT_KPI*r)*Rf[t*R+(r+R/2)];
00184      for(r=0; r<R/2; r++)
00185      fft[r+R/2] = cexp(I*NFFT_KPI*r)*Rf[t*R+r];
00186      */
00187 
00188     for (r = 0; r < S; r++)
00189       fft[r] = Rf[t * S + r] + _Complex_I * NFFT_K(0.0);
00190 
00191     NFFT(fftshift_complex_int)(fft, 1, &S);
00192     FFTW(execute)(my_fftw_plan);
00193     NFFT(fftshift_complex_int)(fft, 1, &S);
00194 
00195     my_infft_plan.y[t * S] = NFFT_K(0.0);
00196     for (r = -S / 2 + 1; r < S / 2; r++)
00197       my_infft_plan.y[t * S + (r + S / 2)] = fft[r + S / 2] / KERNEL(r);
00198   }
00199 
00201   for (k = 0; k < my_nfft_plan.N_total; k++)
00202     my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
00203 
00205   SOLVER(before_loop_complex)(&my_infft_plan);
00206 
00207   if (max_i < 1)
00208   {
00209     l = 1;
00210     for (k = 0; k < my_nfft_plan.N_total; k++)
00211       my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
00212   }
00213   else
00214   {
00215     for (l = 1; l <= max_i; l++)
00216     {
00217       SOLVER(loop_one_step_complex)(&my_infft_plan);
00218       /*if (sqrt(my_infft_plan.dot_r_iter)<=1e-12) break;*/
00219     }
00220   }
00221   /*printf("after %d iteration(s): weighted 2-norm of original residual vector = %g\n",l-1,sqrt(my_infft_plan.dot_r_iter));*/
00222 
00224   for (k = 0; k < my_nfft_plan.N_total; k++)
00225     f[k] = NFFT_M(creal)(my_infft_plan.f_hat_iter[k]);
00226 
00228   FFTW(destroy_plan)(my_fftw_plan);
00229   NFFT(free)(fft);
00230   SOLVER(finalize_complex)(&my_infft_plan);
00231   NFFT(finalize)(&my_nfft_plan);
00232   NFFT(free)(x);
00233   NFFT(free)(w);
00234   return 0;
00235 }
00236 
00239 int main(int argc, char **argv)
00240 {
00241   int (*gridfcn)(); 
00242   int T, S; 
00243   FILE *fp;
00244   int N; 
00245   NFFT_R *Rf, *iRf;
00246   int max_i; 
00248   if (argc != 6)
00249   {
00250     printf("inverse_radon gridfcn N T R max_i\n");
00251     printf("\n");
00252     printf("gridfcn    \"polar\" or \"linogram\" \n");
00253     printf("N          image size NxN            \n");
00254     printf("T          number of slopes          \n");
00255     printf("R          number of offsets         \n");
00256     printf("max_i      number of iterations      \n");
00257     exit(EXIT_FAILURE);
00258   }
00259 
00260   if (strcmp(argv[1], "polar") == 0)
00261     gridfcn = polar_grid;
00262   else
00263     gridfcn = linogram_grid;
00264 
00265   N = atoi(argv[2]);
00266   T = atoi(argv[3]);
00267   S = atoi(argv[4]);
00268   /*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
00269   max_i = atoi(argv[5]);
00270 
00271   Rf = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
00272   iRf = (NFFT_R *) NFFT(malloc)((size_t)(N * N) * (sizeof(NFFT_R)));
00273 
00275   fp = fopen("sinogram_data.bin", "rb");
00276   if (fp == NULL)
00277     return EXIT_FAILURE;
00278   fread(Rf, sizeof(NFFT_R), (size_t)(T * S), fp);
00279   fclose(fp);
00280 
00282   inverse_radon_trafo(gridfcn, T, S, Rf, N, iRf, max_i);
00283 
00285   fp = fopen("output_data.bin", "wb+");
00286   if (fp == NULL)
00287     return EXIT_FAILURE;
00288   fwrite(iRf, sizeof(NFFT_R), (size_t)(N * N), fp);
00289   fclose(fp);
00290 
00292   NFFT(free)(Rf);
00293   NFFT(free)(iRf);
00294 
00295   return EXIT_SUCCESS;
00296 }