NFFT  3.3.0
radon.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 
41 #include <stdio.h>
42 #include <math.h>
43 #include <stdlib.h>
44 #include <string.h>
45 #include <complex.h>
46 
47 #define NFFT_PRECISION_DOUBLE
48 
49 #include "nfft3mp.h"
50 
52 /*#define KERNEL(r) 1.0 */
53 #define KERNEL(r) (NFFT_K(1.0)-NFFT_M(fabs)((NFFT_R)(r))/((NFFT_R)S/2))
54 
58 static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
59 {
60  int t, r;
61  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));
62 
63  for (t = -T / 2; t < T / 2; t++)
64  {
65  for (r = -S / 2; r < S / 2; r++)
66  {
67  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));
68  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));
69  if (r == 0)
70  w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
71  else
72  w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
73  }
74  }
75 
76  return 0;
77 }
78 
82 static int linogram_grid(int T, int S, NFFT_R *x, NFFT_R *w)
83 {
84  int t, r;
85  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));
86 
87  for (t = -T / 2; t < T / 2; t++)
88  {
89  for (r = -S / 2; r < S / 2; r++)
90  {
91  if (t < 0)
92  {
93  x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S);
94  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)
95  / (NFFT_R)(S);
96  }
97  else
98  {
99  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)
100  * (NFFT_R)(r) / (NFFT_R)(S);
101  x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S);
102  }
103  if (r == 0)
104  w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
105  else
106  w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
107  }
108  }
109 
110  return 0;
111 }
112 
116 static int Radon_trafo(int (*gridfcn)(), int T, int S, NFFT_R *f, int NN, NFFT_R *Rf)
117 {
118  int j, k;
119  NFFT(plan) my_nfft_plan;
121  NFFT_C *fft;
122  FFTW(plan) my_fftw_plan;
124  int t, r;
125  NFFT_R *x, *w;
127  int N[2], n[2];
128  int M = T * S;
129 
130  N[0] = NN;
131  n[0] = 2 * N[0];
132  N[1] = NN;
133  n[1] = 2 * N[1];
134 
135  fft = (NFFT_C *) NFFT(malloc)((size_t)(S) * sizeof(NFFT_C));
136  my_fftw_plan = FFTW(plan_dft_1d)(S, fft, fft, FFTW_BACKWARD, FFTW_MEASURE);
137 
138  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
139  if (x == NULL)
140  return EXIT_FAILURE;
141 
142  w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
143  if (w == NULL)
144  return EXIT_FAILURE;
145 
147  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, 4,
148  PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
149  | FFT_OUT_OF_PLACE,
150  FFTW_MEASURE | FFTW_DESTROY_INPUT);
151 
153  gridfcn(T, S, x, w);
154  for (j = 0; j < my_nfft_plan.M_total; j++)
155  {
156  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
157  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
158  }
159 
161  if (my_nfft_plan.flags & PRE_LIN_PSI)
162  NFFT(precompute_lin_psi)(&my_nfft_plan);
163 
164  if (my_nfft_plan.flags & PRE_PSI)
165  NFFT(precompute_psi)(&my_nfft_plan);
166 
167  if (my_nfft_plan.flags & PRE_FULL_PSI)
168  NFFT(precompute_full_psi)(&my_nfft_plan);
169 
171  for (k = 0; k < my_nfft_plan.N_total; k++)
172  my_nfft_plan.f_hat[k] = f[k] + _Complex_I * NFFT_K(0.0);
173 
175  NFFT(trafo)(&my_nfft_plan);
176 
178  for (t = 0; t < T; t++)
179  {
180  fft[0] = NFFT_K(0.0);
181  for (r = -S / 2 + 1; r < S / 2; r++)
182  fft[r + S / 2] = KERNEL(r) * my_nfft_plan.f[t * S + (r + S / 2)];
183 
184  NFFT(fftshift_complex_int)(fft, 1, &S);
185  FFTW(execute)(my_fftw_plan);
186  NFFT(fftshift_complex_int)(fft, 1, &S);
187 
188  for (r = 0; r < S; r++)
189  Rf[t * S + r] = NFFT_M(creal)(fft[r]) / (NFFT_R)(S);
190 
191  /* for(r=0; r<R/2; r++)
192  Rf[t*R+(r+R/2)] = creal(cexp(-I*NFFT_KPI*r)*fft[r]);
193  for(r=0; r<R/2; r++)
194  Rf[t*R+r] = creal(cexp(-I*NFFT_KPI*r)*fft[r+R/2]);
195  */
196  }
197 
199  FFTW(destroy_plan)(my_fftw_plan);
200  NFFT(free)(fft);
201  NFFT(finalize)(&my_nfft_plan);
202  NFFT(free)(x);
203  NFFT(free)(w);
204  return 0;
205 }
206 
209 int main(int argc, char **argv)
210 {
211  int (*gridfcn)();
212  int T, S;
213  FILE *fp;
214  int N;
215  NFFT_R *f, *Rf;
216 
217  if (argc != 5)
218  {
219  printf("radon gridfcn N T R\n");
220  printf("\n");
221  printf("gridfcn \"polar\" or \"linogram\" \n");
222  printf("N image size NxN \n");
223  printf("T number of slopes \n");
224  printf("R number of offsets \n");
225  exit(EXIT_FAILURE);
226  }
227 
228  if (strcmp(argv[1], "polar") == 0)
229  gridfcn = polar_grid;
230  else
231  gridfcn = linogram_grid;
232 
233  N = atoi(argv[2]);
234  T = atoi(argv[3]);
235  S = atoi(argv[4]);
236  /*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
237 
238  f = (NFFT_R *) NFFT(malloc)((size_t)(N * N) * (sizeof(NFFT_R)));
239  Rf = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
240 
242  fp = fopen("input_data.bin", "rb");
243  if (fp == NULL)
244  return EXIT_FAILURE;
245  fread(f, sizeof(NFFT_R), (size_t)(N * N), fp);
246  fclose(fp);
247 
249  Radon_trafo(gridfcn, T, S, f, N, Rf);
250 
252  fp = fopen("sinogram_data.bin", "wb+");
253  if (fp == NULL)
254  return EXIT_FAILURE;
255  fwrite(Rf, sizeof(NFFT_R), (size_t)(T * S), fp);
256  fclose(fp);
257 
259  NFFT(free)(f);
260  NFFT(free)(Rf);
261 
262  return EXIT_SUCCESS;
263 }
static int Radon_trafo(int(*gridfcn)(), int T, int S, NFFT_R *f, int NN, NFFT_R *Rf)
computes the NFFT-based discrete Radon transform of f on the grid given by gridfcn() with T angles an...
Definition: radon.c:116
static int linogram_grid(int T, int S, NFFT_R *x, NFFT_R *w)
generates the points x with weights w for the linogram grid with T slopes and R offsets ...
Definition: radon.c:82
static void fft(int N, int M, int Z, fftw_complex *mem)
fft makes an 1D-ftt for every knot through all layers
static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
generates the points x with weights w for the polar grid with T angles and R offsets ...
Definition: radon.c:58
#define KERNEL(r)
define weights of kernel function for discrete Radon transform
Definition: radon.c:53