NFFT  3.5.3
inverse_radon.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2017 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 
38 #include <stdio.h>
39 #include <math.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <complex.h>
43 
44 #define NFFT_PRECISION_DOUBLE
45 
46 #include "nfft3mp.h"
47 
49 /*#define KERNEL(r) 1.0 */
50 #define KERNEL(r) (NFFT_K(1.0)-NFFT_M(fabs)((NFFT_R)(r))/((NFFT_R)S/2))
51 
55 static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w)
56 {
57  int t, r;
58  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));
59 
60  for (t = -T / 2; t < T / 2; t++)
61  {
62  for (r = -S / 2; r < S / 2; r++)
63  {
64  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));
65  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));
66  if (r == 0)
67  w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
68  else
69  w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
70  }
71  }
72 
73  return 0;
74 }
75 
79 static int linogram_grid(int T, int S, NFFT_R *x, NFFT_R *w)
80 {
81  int t, r;
82  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));
83 
84  for (t = -T / 2; t < T / 2; t++)
85  {
86  for (r = -S / 2; r < S / 2; r++)
87  {
88  if (t < 0)
89  {
90  x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S);
91  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)
92  / (NFFT_R)(S);
93  }
94  else
95  {
96  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)
97  * (NFFT_R)(r) / (NFFT_R)(S);
98  x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S);
99  }
100  if (r == 0)
101  w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
102  else
103  w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
104  }
105  }
106 
107  return 0;
108 }
109 
114 static int inverse_radon_trafo(int (*gridfcn)(), int T, int S, NFFT_R *Rf, int NN, NFFT_R *f,
115  int max_i)
116 {
117  int j, k;
118  NFFT(plan) my_nfft_plan;
119  SOLVER(plan_complex) my_infft_plan;
121  NFFT_C *fft;
122  FFTW(plan) my_fftw_plan;
124  int t, r;
125  NFFT_R *x, *w;
126  int l;
128  int N[2], n[2];
129  int M = T * S;
130 
131  N[0] = NN;
132  n[0] = 2 * N[0];
133  N[1] = NN;
134  n[1] = 2 * N[1];
135 
136  fft = (NFFT_C *) NFFT(malloc)((size_t)(S) * sizeof(NFFT_C));
137  my_fftw_plan = FFTW(plan_dft_1d)(S, fft, fft, FFTW_FORWARD, FFTW_MEASURE);
138 
139  x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R)));
140  if (x == NULL)
141  return EXIT_FAILURE;
142 
143  w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
144  if (w == NULL)
145  return EXIT_FAILURE;
146 
148  NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, 4,
150  FFTW_MEASURE);
151 
153  SOLVER(init_advanced_complex)(&my_infft_plan,
154  (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
155 
157  gridfcn(T, S, x, w);
158  for (j = 0; j < my_nfft_plan.M_total; j++)
159  {
160  my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
161  my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
162  if (j % S)
163  my_infft_plan.w[j] = w[j];
164  else
165  my_infft_plan.w[j] = NFFT_K(0.0);
166  }
167 
169  if (my_nfft_plan.flags & PRE_LIN_PSI)
170  NFFT(precompute_lin_psi)(&my_nfft_plan);
171 
172  if (my_nfft_plan.flags & PRE_PSI)
173  NFFT(precompute_psi)(&my_nfft_plan);
174 
175  if (my_nfft_plan.flags & PRE_FULL_PSI)
176  NFFT(precompute_full_psi)(&my_nfft_plan);
177 
179  for (t = 0; t < T; t++)
180  {
181  /* for(r=0; r<R/2; r++)
182  fft[r] = cexp(I*NFFT_KPI*r)*Rf[t*R+(r+R/2)];
183  for(r=0; r<R/2; r++)
184  fft[r+R/2] = cexp(I*NFFT_KPI*r)*Rf[t*R+r];
185  */
186 
187  for (r = 0; r < S; r++)
188  fft[r] = Rf[t * S + r] + _Complex_I * NFFT_K(0.0);
189 
190  NFFT(fftshift_complex_int)(fft, 1, &S);
191  FFTW(execute)(my_fftw_plan);
192  NFFT(fftshift_complex_int)(fft, 1, &S);
193 
194  my_infft_plan.y[t * S] = NFFT_K(0.0);
195  for (r = -S / 2 + 1; r < S / 2; r++)
196  my_infft_plan.y[t * S + (r + S / 2)] = fft[r + S / 2] / KERNEL(r);
197  }
198 
200  for (k = 0; k < my_nfft_plan.N_total; k++)
201  my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
202 
204  SOLVER(before_loop_complex)(&my_infft_plan);
205 
206  if (max_i < 1)
207  {
208  l = 1;
209  for (k = 0; k < my_nfft_plan.N_total; k++)
210  my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
211  }
212  else
213  {
214  for (l = 1; l <= max_i; l++)
215  {
216  SOLVER(loop_one_step_complex)(&my_infft_plan);
217  /*if (sqrt(my_infft_plan.dot_r_iter)<=1e-12) break;*/
218  }
219  }
220  /*printf("after %d iteration(s): weighted 2-norm of original residual vector = %g\n",l-1,sqrt(my_infft_plan.dot_r_iter));*/
221 
223  for (k = 0; k < my_nfft_plan.N_total; k++)
224  f[k] = NFFT_M(creal)(my_infft_plan.f_hat_iter[k]);
225 
227  FFTW(destroy_plan)(my_fftw_plan);
228  NFFT(free)(fft);
229  SOLVER(finalize_complex)(&my_infft_plan);
230  NFFT(finalize)(&my_nfft_plan);
231  NFFT(free)(x);
232  NFFT(free)(w);
233  return 0;
234 }
235 
238 int main(int argc, char **argv)
239 {
240  int (*gridfcn)();
241  int T, S;
242  FILE *fp;
243  int N;
244  NFFT_R *Rf, *iRf;
245  int max_i;
247  if (argc != 6)
248  {
249  printf("inverse_radon gridfcn N T R max_i\n");
250  printf("\n");
251  printf("gridfcn \"polar\" or \"linogram\" \n");
252  printf("N image size NxN \n");
253  printf("T number of slopes \n");
254  printf("R number of offsets \n");
255  printf("max_i number of iterations \n");
256  exit(EXIT_FAILURE);
257  }
258 
259  if (strcmp(argv[1], "polar") == 0)
260  gridfcn = polar_grid;
261  else
262  gridfcn = linogram_grid;
263 
264  N = atoi(argv[2]);
265  T = atoi(argv[3]);
266  S = atoi(argv[4]);
267  /*printf("N=%d, %s grid with T=%d, R=%d. \n",N,argv[1],T,R);*/
268  max_i = atoi(argv[5]);
269 
270  Rf = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R)));
271  iRf = (NFFT_R *) NFFT(malloc)((size_t)(N * N) * (sizeof(NFFT_R)));
272 
274  fp = fopen("sinogram_data.bin", "rb");
275  if (fp == NULL)
276  return EXIT_FAILURE;
277  fread(Rf, sizeof(NFFT_R), (size_t)(T * S), fp);
278  fclose(fp);
279 
281  inverse_radon_trafo(gridfcn, T, S, Rf, N, iRf, max_i);
282 
284  fp = fopen("output_data.bin", "wb+");
285  if (fp == NULL)
286  return EXIT_FAILURE;
287  fwrite(iRf, sizeof(NFFT_R), (size_t)(N * N), fp);
288  fclose(fp);
289 
291  NFFT(free)(Rf);
292  NFFT(free)(iRf);
293 
294  return EXIT_SUCCESS;
295 }
static int max_i(int a, int b)
max
Definition: fastsum.c:46
static void fft(int N, int M, int Z, fftw_complex *mem)
fft makes an 1D-ftt for every knot through all layers
#define MALLOC_F_HAT
Definition: nfft3.h:188
#define MALLOC_X
Definition: nfft3.h:187
#define PRE_FULL_PSI
Definition: nfft3.h:186
#define PRE_PSI
Definition: nfft3.h:185
#define MALLOC_F
Definition: nfft3.h:189
#define PRE_LIN_PSI
Definition: nfft3.h:183
#define FFTW_INIT
Definition: nfft3.h:191
#define PRE_PHI_HUT
Definition: nfft3.h:181
#define CGNR
Definition: nfft3.h:775
#define PRECOMPUTE_WEIGHT
Definition: nfft3.h:778
#define KERNEL(r)
define weights of kernel function for discrete Radon transform
Definition: inverse_radon.c:50
static int inverse_radon_trafo(int(*gridfcn)(), int T, int S, NFFT_R *Rf, int NN, NFFT_R *f, int max_i)
computes the inverse discrete Radon transform of Rf on the grid given by gridfcn() with T angles and ...
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: inverse_radon.c:55
int main(int argc, char **argv)
simple test program for the inverse discrete Radon transform
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: inverse_radon.c:79