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 19/** 20 * \file polarFFT/polar_fft_test.c 21 * \brief NFFT-based polar FFT and inverse. 22 * 23 * Computes the NFFT-based polar FFT and its inverse for various parameters. 24 * \author Markus Fenn 25 * \date 2006 26 */ 27 28#include <math.h> 29#include <stdlib.h> 30#include <complex.h> 31 32#define @NFFT_PRECISION_MACRO@ 33 34#include "nfft3mp.h" 35 36/** 37 * \defgroup applications_polarFFT_polar polar_fft_test 38 * \ingroup applications_polarFFT 39 * \{ 40 */ 41 42/** Generates the points \f$x_{t,j}\f$ with weights \f$w_{t,j}\f$ 43 * for the polar grid with \f$T\f$ angles and \f$R\f$ offsets. 44 * 45 * The nodes of the polar grid lie on concentric circles around the origin. 46 * They are given for \f$(j,t)^{\top}\in I_R\times I_T\f$ by 47 * a signed radius \f$r_j := \frac{j}{R} \in [-\frac{1}{2},\frac{1}{2})\f$ and 48 * an angle \f$\theta_t := \frac{\pi t}{T} \in [-\frac{\pi}{2},\frac{\pi}{2})\f$ 49 * as 50 * \f[ 51 * x_{t,j} := r_j\left(\cos\theta_t, \sin\theta_t\right)^{\top}\,. 52 * \f] 53 * The total number of nodes is \f$M=TR\f$, 54 * whereas the origin is included multiple times. 55 * 56 * Weights are introduced to compensate for local sampling density variations. 57 * For every point in the sampling set, we associate a small surrounding area. 58 * In case of the polar grid, we choose small ring segments. 59 * The area of such a ring segment around \f$x_{t,j}\f$ (\f$j \ne 0\f$) is 60 * \f[ 61 * w_{t,j} 62 * = \frac{\pi}{2TR^2}\left(\left(|j|+\frac{1}{2}\right)^2- 63 * \left(|j|-\frac{1}{2}\right)^2\right) 64 * = \frac{\pi |j| }{TR^2}\, . 65 * \f] 66 * The area of the small circle of radius \f$\frac{1}{2R}\f$ around the origin is 67 * \f$\frac{\pi}{4R^2}\f$. 68 * Divided by the multiplicity of the origin in the sampling set, we get 69 * \f$w_{t,0} := \frac{\pi}{4TR^2}\f$. 70 * Thus, the sum of all weights is \f$\frac{\pi}{4}(1+\frac{1}{R^2})\f$ and 71 * we divide by this value for normalization. 72 */ 73static int polar_grid(int T, int S, NFFT_R *x, NFFT_R *w) 74{ 75 int t, r; 76 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)); 77 78 for (t = -T / 2; t < T / 2; t++) 79 { 80 for (r = -S / 2; r < S / 2; r++) 81 { 82 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)); 83 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)); 84 if (r == 0) 85 w[(t + T / 2) * S + (r + S / 2)] = NFFT_K(1.0) / NFFT_K(4.0) / W; 86 else 87 w[(t + T / 2) * S + (r + S / 2)] = NFFT_M(fabs)((NFFT_R) r) / W; 88 } 89 } 90 91 return T * S; /** return the number of knots */ 92} 93 94/** discrete polar FFT */ 95static int polar_dft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, 96 int m) 97{ 98 int j, k; /**< index for nodes and frequencies */ 99 NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */ 100 101 NFFT_R *x, *w; /**< knots and associated weights */ 102 103 int N[2], n[2]; 104 int M = T * S; /**< number of knots */ 105 106 N[0] = NN; 107 n[0] = 2 * N[0]; /**< oversampling factor sigma=2 */ 108 N[1] = NN; 109 n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */ 110 111 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R))); 112 if (x == NULL) 113 return EXIT_FAILURE; 114 115 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R))); 116 if (w == NULL) 117 return EXIT_FAILURE; 118 119 /** init two dimensional NFFT plan */ 120 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m, 121 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT, 122 FFTW_MEASURE); 123 124 /** init nodes from polar grid*/ 125 polar_grid(T, S, x, w); 126 for (j = 0; j < my_nfft_plan.M_total; j++) 127 { 128 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0]; 129 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1]; 130 } 131 132 /** init Fourier coefficients from given image */ 133 for (k = 0; k < my_nfft_plan.N_total; k++) 134 my_nfft_plan.f_hat[k] = f_hat[k]; 135 136 /** NDFT-2D */ 137 NFFT(trafo_direct)(&my_nfft_plan); 138 139 /** copy result */ 140 for (j = 0; j < my_nfft_plan.M_total; j++) 141 f[j] = my_nfft_plan.f[j]; 142 143 /** finalise the plans and free the variables */ 144 NFFT(finalize)(&my_nfft_plan); 145 NFFT(free)(x); 146 NFFT(free)(w); 147 148 return EXIT_SUCCESS; 149} 150 151/** NFFT-based polar FFT */ 152static int polar_fft(NFFT_C *f_hat, int NN, NFFT_C *f, int T, int S, 153 int m) 154{ 155 int j, k; /**< index for nodes and freqencies */ 156 NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */ 157 158 NFFT_R *x, *w; /**< knots and associated weights */ 159 160 int N[2], n[2]; 161 int M = T * S; /**< number of knots */ 162 163 N[0] = NN; 164 n[0] = 2 * N[0]; /**< oversampling factor sigma=2 */ 165 N[1] = NN; 166 n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */ 167 168 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R))); 169 if (x == NULL) 170 return EXIT_FAILURE; 171 172 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R))); 173 if (w == NULL) 174 return EXIT_FAILURE; 175 176 /** init two dimensional NFFT plan */ 177 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m, 178 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT 179 | FFT_OUT_OF_PLACE, 180 FFTW_MEASURE | FFTW_DESTROY_INPUT); 181 182 /** init nodes from polar grid*/ 183 polar_grid(T, S, x, w); 184 for (j = 0; j < my_nfft_plan.M_total; j++) 185 { 186 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0]; 187 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1]; 188 } 189 190 /** precompute psi, the entries of the matrix B */ 191 if (my_nfft_plan.flags & PRE_LIN_PSI) 192 NFFT(precompute_lin_psi)(&my_nfft_plan); 193 194 if (my_nfft_plan.flags & PRE_PSI) 195 NFFT(precompute_psi)(&my_nfft_plan); 196 197 if (my_nfft_plan.flags & PRE_FULL_PSI) 198 NFFT(precompute_full_psi)(&my_nfft_plan); 199 200 /** init Fourier coefficients from given image */ 201 for (k = 0; k < my_nfft_plan.N_total; k++) 202 my_nfft_plan.f_hat[k] = f_hat[k]; 203 204 /** NFFT-2D */ 205 NFFT(trafo)(&my_nfft_plan); 206 207 /** copy result */ 208 for (j = 0; j < my_nfft_plan.M_total; j++) 209 f[j] = my_nfft_plan.f[j]; 210 211 /** finalise the plans and free the variables */ 212 NFFT(finalize)(&my_nfft_plan); 213 NFFT(free)(x); 214 NFFT(free)(w); 215 216 return EXIT_SUCCESS; 217} 218 219/** inverse NFFT-based polar FFT */ 220static int inverse_polar_fft(NFFT_C *f, int T, int S, NFFT_C *f_hat, 221 int NN, int max_i, int m) 222{ 223 int j, k; /**< index for nodes and freqencies */ 224 NFFT(plan) my_nfft_plan; /**< plan for the nfft-2D */ 225 SOLVER(plan_complex) my_infft_plan; /**< plan for the inverse nfft */ 226 227 NFFT_R *x, *w; /**< knots and associated weights */ 228 int l; /**< index for iterations */ 229 230 int N[2], n[2]; 231 int M = T * S; /**< number of knots */ 232 233 N[0] = NN; 234 n[0] = 2 * N[0]; /**< oversampling factor sigma=2 */ 235 N[1] = NN; 236 n[1] = 2 * N[1]; /**< oversampling factor sigma=2 */ 237 238 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * T * S) * (sizeof(NFFT_R))); 239 if (x == NULL) 240 return EXIT_FAILURE; 241 242 w = (NFFT_R *) NFFT(malloc)((size_t)(T * S) * (sizeof(NFFT_R))); 243 if (w == NULL) 244 return EXIT_FAILURE; 245 246 /** init two dimensional NFFT plan */ 247 NFFT(init_guru)(&my_nfft_plan, 2, N, M, n, m, 248 PRE_PHI_HUT | PRE_PSI | MALLOC_X | MALLOC_F_HAT | MALLOC_F | FFTW_INIT 249 | FFT_OUT_OF_PLACE, 250 FFTW_MEASURE | FFTW_DESTROY_INPUT); 251 252 /** init two dimensional infft plan */ 253 SOLVER(init_advanced_complex)(&my_infft_plan, 254 (NFFT(mv_plan_complex)*) (&my_nfft_plan), CGNR | PRECOMPUTE_WEIGHT); 255 256 /** init nodes, given samples and weights */ 257 polar_grid(T, S, x, w); 258 for (j = 0; j < my_nfft_plan.M_total; j++) 259 { 260 my_nfft_plan.x[2 * j + 0] = x[2 * j + 0]; 261 my_nfft_plan.x[2 * j + 1] = x[2 * j + 1]; 262 my_infft_plan.y[j] = f[j]; 263 my_infft_plan.w[j] = w[j]; 264 } 265 266 /** precompute psi, the entries of the matrix B */ 267 if (my_nfft_plan.flags & PRE_LIN_PSI) 268 NFFT(precompute_lin_psi)(&my_nfft_plan); 269 270 if (my_nfft_plan.flags & PRE_PSI) 271 NFFT(precompute_psi)(&my_nfft_plan); 272 273 if (my_nfft_plan.flags & PRE_FULL_PSI) 274 NFFT(precompute_full_psi)(&my_nfft_plan); 275 276 /** initialise damping factors */ 277 if (my_infft_plan.flags & PRECOMPUTE_DAMP) 278 for (j = 0; j < my_nfft_plan.N[0]; j++) 279 for (k = 0; k < my_nfft_plan.N[1]; k++) 280 { 281 my_infft_plan.w_hat[j * my_nfft_plan.N[1] + k] = ( 282 NFFT_M(sqrt)( 283 NFFT_M(pow)((NFFT_R) (j - my_nfft_plan.N[0] / 2), NFFT_K(2.0)) 284 + NFFT_M(pow)((NFFT_R) (k - my_nfft_plan.N[1] / 2), NFFT_K(2.0))) 285 > ((NFFT_R) (my_nfft_plan.N[0] / 2)) ? 0 : 1); 286 } 287 288 /** initialise some guess f_hat_0 */ 289 for (k = 0; k < my_nfft_plan.N_total; k++) 290 my_infft_plan.f_hat_iter[k] = NFFT_K(0.0) + _Complex_I * NFFT_K(0.0); 291 292 /** solve the system */ 293 SOLVER(before_loop_complex)(&my_infft_plan); 294 295 if (max_i < 1) 296 { 297 l = 1; 298 for (k = 0; k < my_nfft_plan.N_total; k++) 299 my_infft_plan.f_hat_iter[k] = my_infft_plan.p_hat_iter[k]; 300 } 301 else 302 { 303 for (l = 1; l <= max_i; l++) 304 { 305 SOLVER(loop_one_step_complex)(&my_infft_plan); 306 } 307 } 308 309 /** copy result */ 310 for (k = 0; k < my_nfft_plan.N_total; k++) 311 f_hat[k] = my_infft_plan.f_hat_iter[k]; 312 313 /** finalise the plans and free the variables */ 314 SOLVER(finalize_complex)(&my_infft_plan); 315 NFFT(finalize)(&my_nfft_plan); 316 NFFT(free)(x); 317 NFFT(free)(w); 318 319 return EXIT_SUCCESS; 320} 321 322/** test program for various parameters */ 323int main(int argc, char **argv) 324{ 325 int N; /**< mpolar FFT size NxN */ 326 int T, S; /**< number of directions/offsets */ 327 int M; /**< number of knots of mpolar grid */ 328 NFFT_R *x, *w; /**< knots and associated weights */ 329 NFFT_C *f_hat, *f, *f_direct, *f_tilde; 330 int k; 331 int max_i; /**< number of iterations */ 332 int m = 1; 333 NFFT_R temp1, temp2, E_max = NFFT_K(0.0); 334 FILE *fp1, *fp2; 335 char filename[30]; 336 337 if (argc != 4) 338 { 339 printf("polar_fft_test N T R \n"); 340 printf("\n"); 341 printf("N polar FFT of size NxN \n"); 342 printf("T number of slopes \n"); 343 printf("R number of offsets \n"); 344 exit(EXIT_FAILURE); 345 } 346 347 N = atoi(argv[1]); 348 T = atoi(argv[2]); 349 S = atoi(argv[3]); 350 printf("N=%d, polar grid with T=%d, R=%d => ", N, T, S); 351 352 x = (NFFT_R *) NFFT(malloc)((size_t)(2 * 5 * (T / 2) * (S / 2)) * (sizeof(NFFT_R))); 353 w = (NFFT_R *) NFFT(malloc)((size_t)(5 * (T / 2) * (S / 2)) * (sizeof(NFFT_R))); 354 355 f_hat = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N)); 356 f = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(T * S)); 357 f_direct = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(T * S)); 358 f_tilde = (NFFT_C *) NFFT(malloc)(sizeof(NFFT_C) * (size_t)(N * N)); 359 360 /** generate knots of mpolar grid */ 361 M = polar_grid(T, S, x, w); 362 printf("M=%d.\n", M); 363 364 /** load data */ 365 fp1 = fopen("input_data_r.dat", "r"); 366 fp2 = fopen("input_data_i.dat", "r"); 367 if (fp1 == NULL) 368 return (-1); 369 for (k = 0; k < N * N; k++) 370 { 371 fscanf(fp1, NFFT__FR__ " ", &temp1); 372 fscanf(fp2, NFFT__FR__ " ", &temp2); 373 f_hat[k] = temp1 + _Complex_I * temp2; 374 } 375 fclose(fp1); 376 fclose(fp2); 377 378 /** direct polar FFT */ 379 polar_dft(f_hat, N, f_direct, T, S, m); 380 // polar_fft(f_hat,N,f_direct,T,R,12); 381 382 /** Test of the polar FFT with different m */ 383 printf("\nTest of the polar FFT: \n"); 384 fp1 = fopen("polar_fft_error.dat", "w+"); 385 for (m = 1; m <= 12; m++) 386 { 387 /** fast polar FFT */ 388 polar_fft(f_hat, N, f, T, S, m); 389 390 /** compute error of fast polar FFT */ 391 E_max = NFFT(error_l_infty_complex)(f_direct, f, M); 392 printf("m=%2d: E_max = %" NFFT__FES__ "\n", m, E_max); 393 fprintf(fp1, "%" NFFT__FES__ "\n", E_max); 394 } 395 fclose(fp1); 396 397 /** Test of the inverse polar FFT for different m in dependece of the iteration number*/ 398 for (m = 3; m <= 9; m += 3) 399 { 400 printf("\nTest of the inverse polar FFT for m=%d: \n", m); 401 sprintf(filename, "polar_ifft_error%d.dat", m); 402 fp1 = fopen(filename, "w+"); 403 for (max_i = 0; max_i <= 100; max_i += 10) 404 { 405 /** inverse polar FFT */ 406 inverse_polar_fft(f_direct, T, S, f_tilde, N, max_i, m); 407 408 /** compute maximum relative error */ 409 /* E_max=0.0; 410 for(k=0;k<N*N;k++) 411 { 412 temp = cabs((f_hat[k]-f_tilde[k])/f_hat[k]); 413 if (temp>E_max) E_max=temp; 414 } 415 */ 416 E_max = NFFT(error_l_infty_complex)(f_hat, f_tilde, N * N); 417 printf("%3d iterations: E_max = %" NFFT__FES__ "\n", max_i, E_max); 418 fprintf(fp1, "%" NFFT__FES__ "\n", E_max); 419 } 420 fclose(fp1); 421 } 422 423 /** free the variables */ 424 NFFT(free)(x); 425 NFFT(free)(w); 426 NFFT(free)(f_hat); 427 NFFT(free)(f); 428 NFFT(free)(f_direct); 429 NFFT(free)(f_tilde); 430 431 return EXIT_SUCCESS; 432} 433/* \} */ 434