NFFT  3.3.1
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 
00039 #include <stdio.h>
00040 #include <math.h>
00041 #include <stdlib.h>
00042 #include <string.h>
00043 #include <complex.h>
00044 
00045 #define NFFT_PRECISION_DOUBLE
00046 
00047 #include "nfft3mp.h"
00048 
00050 /*#define KERNEL(r) 1.0 */
00051 #define KERNEL(r) (NFFT_K(1.0)-NFFT_M(fabs)((NFFT_R)(r))/((NFFT_R)S/2))
00052 
00056 static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
00057 {
00058   int t, r;
00059   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));
00060 
00061   for (t = -T / 2; t < T / 2; t++)
00062   {
00063     for (r = -S / 2; r < S / 2; r++)
00064     {
00065       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));
00066       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));
00067       if (r == 0)
00068         w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
00069       else
00070         w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
00071     }
00072   }
00073 
00074   return 0;
00075 }
00076 
00080 static int linogram_grid(int T, int S, NFFT_R *x, NFFT_R *w)
00081 {
00082   int t, r;
00083   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));
00084 
00085   for (t = -T / 2; t < T / 2; t++)
00086   {
00087     for (r = -S / 2; r < S / 2; r++)
00088     {
00089       if (t < 0)
00090       {
00091         x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S);
00092         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)
00093             / (NFFT_R)(S);
00094       }
00095       else
00096       {
00097         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)
00098             * (NFFT_R)(r) / (NFFT_R)(S);
00099         x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S);
00100       }
00101       if (r == 0)
00102         w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
00103       else
00104         w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
00105     }
00106   }
00107 
00108   return 0;
00109 }
00110 
00114 static int Radon_trafo(int (*gridfcn)(), int T, int S, NFFT_R *f, int NN, NFFT_R *Rf)
00115 {
00116   int j, k; 
00117   NFFT(plan) my_nfft_plan; 
00119   NFFT_C *fft; 
00120   FFTW(plan) my_fftw_plan; 
00122   int t, r; 
00123   NFFT_R *x, *w; 
00125   int N[2], n[2];
00126   int M = T * S;
00127 
00128   N[0] = NN;
00129   n[0] = 2 * N[0];
00130   N[1] = NN;
00131   n[1] = 2 * N[1];
00132 
00133   fft = (NFFT_C *) NFFT(malloc)((size_t)(S) * sizeof(NFFT_C));
00134   my_fftw_plan = FFTW(plan_dft_1d)(S, fft, fft, FFTW_BACKWARD, FFTW_MEASURE);
00135 
00136   x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
00137   if (x == NULL)
00138     return EXIT_FAILURE;
00139 
00140   w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
00141   if (w == NULL)
00142     return EXIT_FAILURE;
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
00147           | FFT_OUT_OF_PLACE,
00148       FFTW_MEASURE | FFTW_DESTROY_INPUT);
00149 
00151   gridfcn(T, S, 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.flags & PRE_LIN_PSI)
00160     NFFT(precompute_lin_psi)(&my_nfft_plan);
00161 
00162   if (my_nfft_plan.flags & PRE_PSI)
00163     NFFT(precompute_psi)(&my_nfft_plan);
00164 
00165   if (my_nfft_plan.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 * NFFT_K(0.0);
00171 
00173   NFFT(trafo)(&my_nfft_plan);
00174 
00176   for (t = 0; t < T; t++)
00177   {
00178     fft[0] = NFFT_K(0.0);
00179     for (r = -S / 2 + 1; r < S / 2; r++)
00180       fft[r + S / 2] = KERNEL(r) * my_nfft_plan.f[t * S + (r + S / 2)];
00181 
00182     NFFT(fftshift_complex_int)(fft, 1, &S);
00183     FFTW(execute)(my_fftw_plan);
00184     NFFT(fftshift_complex_int)(fft, 1, &S);
00185 
00186     for (r = 0; r < S; r++)
00187       Rf[t * S + r] = NFFT_M(creal)(fft[r]) / (NFFT_R)(S);
00188 
00189     /*    for(r=0; r<R/2; r++)
00190      Rf[t*R+(r+R/2)] = creal(cexp(-I*NFFT_KPI*r)*fft[r]);
00191      for(r=0; r<R/2; r++)
00192      Rf[t*R+r] = creal(cexp(-I*NFFT_KPI*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, S; 
00211   FILE *fp;
00212   int N; 
00213   NFFT_R *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(EXIT_FAILURE);
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   S = atoi(argv[4]);
00234   /*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
00235 
00236   f = (NFFT_R *) NFFT(malloc)((size_t)(N * N) * (sizeof(NFFT_R)));
00237   Rf = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
00238 
00240   fp = fopen("input_data.bin", "rb");
00241   if (fp == NULL)
00242     return EXIT_FAILURE;
00243   fread(f, sizeof(NFFT_R), (size_t)(N * N), fp);
00244   fclose(fp);
00245 
00247   Radon_trafo(gridfcn, T, S, f, N, Rf);
00248 
00250   fp = fopen("sinogram_data.bin", "wb+");
00251   if (fp == NULL)
00252     return EXIT_FAILURE;
00253   fwrite(Rf, sizeof(NFFT_R), (size_t)(T * S), fp);
00254   fclose(fp);
00255 
00257   NFFT(free)(f);
00258   NFFT(free)(Rf);
00259 
00260   return EXIT_SUCCESS;
00261 }