46 #define NFFT_PRECISION_DOUBLE
52 #define KERNEL(r) (NFFT_K(1.0)-NFFT_M(fabs)((NFFT_R)(r))/((NFFT_R)S/2))
57 static int polar_grid(
int T,
int S, NFFT_R *x, NFFT_R *w)
60 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 for (t = -T / 2; t < T / 2; t++)
64 for (r = -S / 2; r < S / 2; r++)
66 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));
67 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 w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
71 w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
84 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 for (t = -T / 2; t < T / 2; t++)
88 for (r = -S / 2; r < S / 2; r++)
92 x[2 * ((t + T / 2) * S + (r + S / 2)) + 0] = (NFFT_R) r / (NFFT_R)(S);
93 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)
98 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)
99 * (NFFT_R)(r) / (NFFT_R)(S);
100 x[2 * ((t + T / 2) * S + (r + S / 2)) + 1] = (NFFT_R) r / (NFFT_R)(S);
103 w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
105 w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
120 NFFT(plan) my_nfft_plan;
121 SOLVER(plan_complex) my_infft_plan;
124 FFTW(plan) my_fftw_plan;
138 fft = (NFFT_C *) NFFT(malloc)((size_t)(S) *
sizeof(NFFT_C));
139 my_fftw_plan = FFTW(plan_dft_1d)(S,
fft,
fft, FFTW_FORWARD, FFTW_MEASURE);
141 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (
sizeof(NFFT_R)));
145 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (
sizeof(NFFT_R)));
150 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, 4,
151 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
153 FFTW_MEASURE | FFTW_DESTROY_INPUT);
156 SOLVER(init_advanced_complex)(&my_infft_plan,
157 (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
161 for (j = 0; j < my_nfft_plan.M_total; j++)
163 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
164 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
166 my_infft_plan.w[j] = w[j];
168 my_infft_plan.w[j] = NFFT_K(0.0);
172 if (my_nfft_plan.flags & PRE_LIN_PSI)
173 NFFT(precompute_lin_psi)(&my_nfft_plan);
175 if (my_nfft_plan.flags & PRE_PSI)
176 NFFT(precompute_psi)(&my_nfft_plan);
178 if (my_nfft_plan.flags & PRE_FULL_PSI)
179 NFFT(precompute_full_psi)(&my_nfft_plan);
182 for (t = 0; t < T; t++)
190 for (r = 0; r < S; r++)
191 fft[r] = Rf[t * S + r] + _Complex_I * NFFT_K(0.0);
193 NFFT(fftshift_complex_int)(
fft, 1, &S);
194 FFTW(execute)(my_fftw_plan);
195 NFFT(fftshift_complex_int)(
fft, 1, &S);
197 my_infft_plan.y[t * S] = NFFT_K(0.0);
198 for (r = -S / 2 + 1; r < S / 2; r++)
199 my_infft_plan.y[t * S + (r + S / 2)] = fft[r + S / 2] /
KERNEL(r);
203 for (k = 0; k < my_nfft_plan.N_total; k++)
204 my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
207 SOLVER(before_loop_complex)(&my_infft_plan);
212 for (k = 0; k < my_nfft_plan.N_total; k++)
213 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
217 for (l = 1; l <=
max_i; l++)
219 SOLVER(loop_one_step_complex)(&my_infft_plan);
226 for (k = 0; k < my_nfft_plan.N_total; k++)
227 f[k] = NFFT_M(creal)(my_infft_plan.f_hat_iter[k]);
230 FFTW(destroy_plan)(my_fftw_plan);
232 SOLVER(finalize_complex)(&my_infft_plan);
233 NFFT(finalize)(&my_nfft_plan);
241 int main(
int argc,
char **argv)
252 printf(
"inverse_radon gridfcn N T R max_i\n");
254 printf(
"gridfcn \"polar\" or \"linogram\" \n");
255 printf(
"N image size NxN \n");
256 printf(
"T number of slopes \n");
257 printf(
"R number of offsets \n");
258 printf(
"max_i number of iterations \n");
262 if (strcmp(argv[1],
"polar") == 0)
271 max_i = atoi(argv[5]);
273 Rf = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (
sizeof(NFFT_R)));
274 iRf = (NFFT_R *) NFFT(malloc)((size_t)(N * N) * (
sizeof(NFFT_R)));
277 fp = fopen(
"sinogram_data.bin",
"rb");
280 fread(Rf,
sizeof(NFFT_R), (
size_t)(T * S), fp);
287 fp = fopen(
"output_data.bin",
"wb+");
290 fwrite(iRf,
sizeof(NFFT_R), (
size_t)(N * N), fp);
static int max_i(int a, int b)
max
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 ...
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 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 ...
#define KERNEL(r)
define weights of kernel function for discrete Radon transform
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 ...