34 #define NFFT_PRECISION_DOUBLE
44 NFFT_R GLOBAL_elapsed_time;
52 NFFT_R W = (NFFT_R) T * (((NFFT_R) rr / NFFT_K(2.0)) * ((NFFT_R) rr / NFFT_K(2.0)) + NFFT_K(1.0) / NFFT_K(4.0));
54 for (t = -T / 2; t < T / 2; t++)
56 for (r = -rr / 2; r < rr / 2; r++)
60 x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 0] = (NFFT_R) (r) / (NFFT_R)(rr);
61 x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 1] = NFFT_K(4.0) * ((NFFT_R)(t) + (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T)
62 * (NFFT_R)(r) / (NFFT_R)(rr);
66 x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 0] = -NFFT_K(4.0) * ((NFFT_R)(t) - (NFFT_R)(T) / NFFT_K(4.0)) / (NFFT_R)(T)
67 * (NFFT_R)(r) / (NFFT_R)(rr);
68 x[2 * ((t + T / 2) * rr + (r + rr / 2)) + 1] = (NFFT_R) r / (NFFT_R)(rr);
71 w[(t + T / 2) * rr + (r + rr / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W;
73 w[(t + T / 2) * rr + (r + rr / 2)] = NFFT_M(fabs)((NFFT_R) r) / W;
81 static int linogram_dft(NFFT_C *f_hat,
int NN, NFFT_C *f,
int T,
int rr,
int m)
85 NFFT(plan) my_nfft_plan;
97 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (
sizeof(NFFT_R)));
101 w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (
sizeof(NFFT_R)));
106 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
107 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
109 FFTW_MEASURE | FFTW_DESTROY_INPUT);
113 for (j = 0; j < my_nfft_plan.M_total; j++)
115 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
116 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
120 for (k = 0; k < my_nfft_plan.N_total; k++)
121 my_nfft_plan.f_hat[k] = f_hat[k];
124 t0 = NFFT(clock_gettime_seconds)();
125 NFFT(trafo_direct)(&my_nfft_plan);
126 t1 = NFFT(clock_gettime_seconds)();
127 GLOBAL_elapsed_time = (t1 - t0);
130 for (j = 0; j < my_nfft_plan.M_total; j++)
131 f[j] = my_nfft_plan.f[j];
134 NFFT(finalize)(&my_nfft_plan);
142 static int linogram_fft(NFFT_C *f_hat,
int NN, NFFT_C *f,
int T,
int rr,
int m)
146 NFFT(plan) my_nfft_plan;
158 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (
sizeof(NFFT_R)));
162 w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (
sizeof(NFFT_R)));
167 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
168 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
170 FFTW_MEASURE | FFTW_DESTROY_INPUT);
174 for (j = 0; j < my_nfft_plan.M_total; j++)
176 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
177 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
181 if (my_nfft_plan.flags & PRE_LIN_PSI)
182 NFFT(precompute_lin_psi)(&my_nfft_plan);
184 if (my_nfft_plan.flags & PRE_PSI)
185 NFFT(precompute_psi)(&my_nfft_plan);
187 if (my_nfft_plan.flags & PRE_FULL_PSI)
188 NFFT(precompute_full_psi)(&my_nfft_plan);
191 for (k = 0; k < my_nfft_plan.N_total; k++)
192 my_nfft_plan.f_hat[k] = f_hat[k];
195 t0 = NFFT(clock_gettime_seconds)();
196 NFFT(trafo)(&my_nfft_plan);
197 t1 = NFFT(clock_gettime_seconds)();
198 GLOBAL_elapsed_time = (t1 - t0);
201 for (j = 0; j < my_nfft_plan.M_total; j++)
202 f[j] = my_nfft_plan.f[j];
205 NFFT(finalize)(&my_nfft_plan);
218 NFFT(plan) my_nfft_plan;
219 SOLVER(plan_complex) my_infft_plan;
232 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * rr) * (
sizeof(NFFT_R)));
236 w = (NFFT_R *) NFFT(malloc)((size_t)(T * rr) * (
sizeof(NFFT_R)));
241 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m,
242 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT
244 FFTW_MEASURE | FFTW_DESTROY_INPUT);
247 SOLVER(init_advanced_complex)(&my_infft_plan,
248 (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT);
252 for (j = 0; j < my_nfft_plan.M_total; j++)
254 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0];
255 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1];
256 my_infft_plan.y[j] = f[j];
257 my_infft_plan.w[j] = w[j];
261 if (my_nfft_plan.flags & PRE_LIN_PSI)
262 NFFT(precompute_lin_psi)(&my_nfft_plan);
264 if (my_nfft_plan.flags & PRE_PSI)
265 NFFT(precompute_psi)(&my_nfft_plan);
267 if (my_nfft_plan.flags & PRE_FULL_PSI)
268 NFFT(precompute_full_psi)(&my_nfft_plan);
271 if (my_infft_plan.flags & PRECOMPUTE_DAMP)
272 for (j = 0; j < my_nfft_plan.N[0]; j++)
273 for (k = 0; k < my_nfft_plan.N[1]; k++)
275 my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = (
277 NFFT_M(pow)((NFFT_R)(j - my_nfft_plan.N[0] / 2), NFFT_K(2.0))
278 + NFFT_M(pow)((NFFT_R)(k - my_nfft_plan.N[1] / 2), NFFT_K(2.0)))
279 > (NFFT_R)(my_nfft_plan.N[0] / 2) ? 0 : 1);
283 for (k = 0; k < my_nfft_plan.N_total; k++)
284 my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0);
286 t0 = NFFT(clock_gettime_seconds)();
288 SOLVER(before_loop_complex)(&my_infft_plan);
293 for (k = 0; k < my_nfft_plan.N_total; k++)
294 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k];
298 for (l = 1; l <=
max_i; l++)
300 SOLVER(loop_one_step_complex)(&my_infft_plan);
304 t1 = NFFT(clock_gettime_seconds)();
305 GLOBAL_elapsed_time = (t1 - t0);
308 for (k = 0; k < my_nfft_plan.N_total; k++)
309 f_hat[k] = my_infft_plan.f_hat_iter[k];
312 SOLVER(finalize_complex)(&my_infft_plan);
313 NFFT(finalize)(&my_nfft_plan);
324 FFTW(plan) my_fftw_plan;
327 NFFT_R t_fft, t_dft_linogram;
329 f_hat = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(N * N));
330 f = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)((T * rr / 4) * 5));
332 my_fftw_plan = FFTW(plan_dft_2d)(N, N, f_hat, f, FFTW_BACKWARD, FFTW_MEASURE);
334 for (k = 0; k < N * N; k++)
335 f_hat[k] = NFFT(drand48)() + _Complex_I * NFFT(drand48)();
337 t0 = NFFT(clock_gettime_seconds)();
338 for (m = 0; m < 65536 / N; m++)
340 FFTW(execute)(my_fftw_plan);
342 f_hat[2] = NFFT_K(2.0) * f_hat[0];
344 t1 = NFFT(clock_gettime_seconds)();
345 GLOBAL_elapsed_time = (t1 - t0);
346 t_fft = (NFFT_R)(N) * GLOBAL_elapsed_time / NFFT_K(65536.0);
351 t_dft_linogram = GLOBAL_elapsed_time;
354 for (m = 3; m <= 9; m += 3)
356 if ((m == 3) && (N < 256))
357 fprintf(fp,
"%d\t&\t&\t%1.1" NFFT__FES__
"&\t%1.1" NFFT__FES__
"&\t%d\t", N, t_fft, t_dft_linogram,
360 fprintf(fp,
"%d\t&\t&\t%1.1" NFFT__FES__
"&\t &\t%d\t", N, t_fft, m);
362 fprintf(fp,
" \t&\t&\t &\t &\t%d\t", m);
364 printf(
"N=%d\tt_fft=%1.1" NFFT__FES__
"\tt_dft_linogram=%1.1" NFFT__FES__
"\tm=%d\t", N, t_fft,
368 fprintf(fp,
"%1.1" NFFT__FES__
"&\t", GLOBAL_elapsed_time);
369 printf(
"t_linogram=%1.1" NFFT__FES__
"\t", GLOBAL_elapsed_time);
372 fprintf(fp,
"%1.1" NFFT__FES__
"\\\\\\hline\n", GLOBAL_elapsed_time);
374 fprintf(fp,
"%1.1" NFFT__FES__
"\\\\\n", GLOBAL_elapsed_time);
375 printf(
"t_ilinogram=%1.1" NFFT__FES__
"\n", GLOBAL_elapsed_time);
387 int main(
int argc,
char **argv)
393 NFFT_C *f_hat, *f, *f_direct, *f_tilde;
397 NFFT_R temp1, temp2, E_max = NFFT_K(0.0);
404 printf(
"linogram_fft_test N T R \n");
406 printf(
"N linogram FFT of size NxN \n");
407 printf(
"T number of slopes \n");
408 printf(
"R number of offsets \n");
411 printf(
"\nHence, comparison FFTW, linogram FFT and inverse linogram FFT\n");
412 fp1 = fopen(
"linogram_comparison_fft.dat",
"w");
415 for (logN = 4; logN <= 8; logN++)
417 3 * (
int)(1U << (logN - 1)));
426 printf(
"N=%d, linogram grid with T=%d, R=%d => ", N, T, rr);
428 x = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * rr / 2) * (
sizeof(NFFT_R)));
429 w = (NFFT_R *) NFFT(malloc)((size_t)(5 * T * rr / 4) * (
sizeof(NFFT_R)));
431 f_hat = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(N * N));
432 f = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(5 * T * rr / 4));
433 f_direct = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(5 * T * rr / 4));
434 f_tilde = (NFFT_C *) NFFT(malloc)(
sizeof(NFFT_C) * (
size_t)(N * N));
438 printf(
"M=%d.\n", M);
441 fp1 = fopen(
"input_data_r.dat",
"r");
442 fp2 = fopen(
"input_data_i.dat",
"r");
443 if ((fp1 == NULL) || (fp2 == NULL))
445 for (k = 0; k < N * N; k++)
447 fscanf(fp1, NFFT__FR__
" ", &temp1);
448 fscanf(fp2, NFFT__FR__
" ", &temp2);
449 f_hat[k] = temp1 + _Complex_I * temp2;
459 printf(
"\nTest of the linogram FFT: \n");
460 fp1 = fopen(
"linogram_fft_error.dat",
"w+");
461 for (m = 1; m <= 12; m++)
467 E_max = NFFT(error_l_infty_complex)(f_direct, f, M);
470 printf(
"m=%2d: E_max = %" NFFT__FES__
"\n", m, E_max);
471 fprintf(fp1,
"%" NFFT__FES__
"\n", E_max);
476 for (m = 3; m <= 9; m += 3)
478 printf(
"\nTest of the inverse linogram FFT for m=%d: \n", m);
479 sprintf(filename,
"linogram_ifft_error%d.dat", m);
480 fp1 = fopen(filename,
"w+");
481 for (max_i = 0; max_i <= 20; max_i += 2)
487 E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N);
488 printf(
"%3d iterations: E_max = %" NFFT__FES__
"\n", max_i, E_max);
489 fprintf(fp1,
"%" NFFT__FES__
"\n", E_max);
499 NFFT(free)(f_direct);
static int max_i(int a, int b)
max
static int inverse_linogram_fft(NFFT_C *f, int T, int rr, NFFT_C *f_hat, int NN, int max_i, int m)
NFFT-based inverse pseudo-polar FFT.
static int linogram_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int rr, int m)
discrete pseudo-polar FFT
static int comparison_fft(FILE *fp, int N, int T, int rr)
Comparison of the FFTW, linogram FFT, and inverse linogram FFT.
static int linogram_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int rr, int m)
NFFT-based pseudo-polar FFT.
static int linogram_grid(int T, int rr, NFFT_R *x, NFFT_R *w)
Generates the points x with weights w for the linogram grid with T slopes and R offsets.