NFFT 3.5.3alpha
nsfft.c
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#include "config.h"
19
20#include <stdio.h>
21#include <math.h>
22#include <string.h>
23#include <stdlib.h>
24#ifdef HAVE_COMPLEX_H
25#include <complex.h>
26#endif
27#include "nfft3.h"
28#include "infft.h"
29
30#define NSFTT_DISABLE_TEST
31
32/* computes a 2d ndft by 1d nfft along the dimension 1 times
33 1d ndft along dimension 0
34*/
35static void short_nfft_trafo_2d(nfft_plan* ths, nfft_plan* plan_1d)
36{
37 int j,k0;
38 double omega;
39
40 for(j=0;j<ths->M_total;j++)
41 {
42 ths->f[j]= 0;
43 plan_1d->x[j] = ths->x[ths->d * j + 1];
44 }
45
46 for(k0=0;k0<ths->N[0];k0++) /* for shorties */
47 {
48 plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
49
50 nfft_trafo(plan_1d);
51
52 for(j=0;j<ths->M_total;j++)
53 {
54 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
55 ths->f[j] += plan_1d->f[j] * cexp( - I*2*KPI*omega);
56 }
57 }
58}
59
60static void short_nfft_adjoint_2d(nfft_plan* ths, nfft_plan* plan_1d)
61{
62 int j,k0;
63 double omega;
64
65 for(j=0;j<ths->M_total;j++)
66 plan_1d->x[j] = ths->x[ths->d * j + 1];
67
68 for(k0=0;k0<ths->N[0];k0++) /* for shorties */
69 {
70 for(j=0;j<ths->M_total;j++)
71 {
72 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
73 plan_1d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
74 }
75
76 plan_1d->f_hat = ths->f_hat + k0*ths->N[1];
77
78 nfft_adjoint(plan_1d);
79 }
80}
81
82/* computes a 3d ndft by 1d nfft along the dimension 2 times
83 2d ndft along dimension 0,1
84*/
85static void short_nfft_trafo_3d_1(nfft_plan* ths, nfft_plan* plan_1d)
86{
87 int j,k0,k1;
88 double omega;
89
90 for(j=0;j<ths->M_total;j++)
91 {
92 ths->f[j] = 0;
93 plan_1d->x[j] = ths->x[ths->d * j + 2];
94 }
95
96 for(k0=0;k0<ths->N[0];k0++) /* for shorties */
97 for(k1=0;k1<ths->N[1];k1++)
98 {
99 plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
100
101 nfft_trafo(plan_1d);
102
103 for(j=0;j<ths->M_total;j++)
104 {
105 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0]
106 + ((double)(k1 - ths->N[1]/2)) * ths->x[ths->d * j + 1];
107 ths->f[j] += plan_1d->f[j] * cexp( - I*2*KPI*omega);
108 }
109 }
110}
111
112static void short_nfft_adjoint_3d_1(nfft_plan* ths, nfft_plan* plan_1d)
113{
114 int j,k0,k1;
115 double omega;
116
117 for(j=0;j<ths->M_total;j++)
118 plan_1d->x[j] = ths->x[ths->d * j + 2];
119
120 for(k0=0;k0<ths->N[0];k0++) /* for shorties */
121 for(k1=0;k1<ths->N[1];k1++)
122 {
123 for(j=0;j<ths->M_total;j++)
124 {
125 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0]
126 + ((double)(k1 - ths->N[1]/2)) * ths->x[ths->d * j + 1];
127 plan_1d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
128 }
129
130 plan_1d->f_hat = ths->f_hat + (k0*ths->N[1]+k1)*ths->N[2];
131
132 nfft_adjoint(plan_1d);
133 }
134}
135
136/* computes a 3d ndft by 2d nfft along the dimension 1,2 times
137 1d ndft along dimension 0
138*/
139static void short_nfft_trafo_3d_2(nfft_plan* ths, nfft_plan* plan_2d)
140{
141 int j,k0;
142 double omega;
143
144 for(j=0;j<ths->M_total;j++)
145 {
146 ths->f[j] = 0;
147 plan_2d->x[2*j+0] = ths->x[ths->d * j + 1];
148 plan_2d->x[2*j+1] = ths->x[ths->d * j + 2];
149 }
150
151 for(k0=0;k0<ths->N[0];k0++) /* for shorties */
152 {
153 plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
154
155 nfft_trafo(plan_2d);
156
157 for(j=0;j<ths->M_total;j++)
158 {
159 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
160 ths->f[j] += plan_2d->f[j] * cexp( - I*2*KPI*omega);
161 }
162 }
163}
164
165static void short_nfft_adjoint_3d_2(nfft_plan* ths, nfft_plan* plan_2d)
166{
167 int j,k0;
168 double omega;
169
170 for(j=0;j<ths->M_total;j++)
171 {
172 plan_2d->x[2*j+0] = ths->x[ths->d * j + 1];
173 plan_2d->x[2*j+1] = ths->x[ths->d * j + 2];
174 }
175
176 for(k0=0;k0<ths->N[0];k0++) /* for shorties */
177 {
178 for(j=0;j<ths->M_total;j++)
179 {
180 omega = ((double)(k0 - ths->N[0]/2)) * ths->x[ths->d * j + 0];
181 plan_2d->f[j] = ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
182 }
183
184 plan_2d->f_hat = ths->f_hat + k0*ths->N[1]*ths->N[2];
185
186 nfft_adjoint(plan_2d);
187 }
188}
189
190/*---------------------------------------------------------------------------*/
191
192#ifdef GAUSSIAN
193static int index_sparse_to_full_direct_2d(int J, int k)
194{
195 int N=X(exp2i)(J+2); /* number of full coeffs */
196 int N_B=X(exp2i)(J); /* number in each sparse block */
197
198 int j=k/N_B; /* consecutive number of Block */
199 int r=j/4; /* level of block */
200
201 int i, o, a, b,s,l,m1,m2;
202 int k1,k2;
203
204 if (k>=(J+4)*X(exp2i)(J+1))
205 {
206 printf("Fehler!\n");
207 return(-1);
208 }
209 else
210 {
211 if (r>(J+1)/2) /* center block */
212 {
213 i=k-4*((J+1)/2+1)*N_B;
214 a=X(exp2i)(J/2+1);
215 m1=i/a;
216 m2=i%a;
217 k1=N/2-a/2+m1;
218 k2=N/2-a/2+m2;
219 }
220 else /* no center block */
221 {
222 i=k-j*N_B; /* index in specific block */
223 o=j%4; /* kind of specific block */
224 a=X(exp2i)(r);
225 b=X(exp2i)(J-r);
226 l=MAX(a,b); /* long dimension of block */
227 s=MIN(a,b); /* short dimension of block */
228 m1=i/l;
229 m2=i%l;
230
231 switch(o)
232 {
233 case 0:
234 {
235 k1=N/2-a/2 ;
236 k2=N/2+ b ;
237
238 if (b>=a)
239 {
240 k1+=m1;
241 k2+=m2;
242 }
243 else
244 {
245 k1+=m2;
246 k2+=m1;
247 }
248 break;
249 }
250 case 1:
251 {
252 k1=N/2+ b ;
253 k2=N/2-a/2 ;
254
255 if (b>a)
256 {
257 k1+=m2;
258 k2+=m1;
259 }
260 else
261 {
262 k1+=m1;
263 k2+=m2;
264 }
265 break;
266 }
267 case 2:
268 {
269 k1=N/2-a/2 ;
270 k2=N/2-2*b ;
271
272 if (b>=a)
273 {
274 k1+=m1;
275 k2+=m2;
276 }
277 else
278 {
279 k1+=m2;
280 k2+=m1;
281 }
282 break;
283 }
284 case 3:
285 {
286 k1=N/2-2*b ;
287 k2=N/2-a/2 ;
288
289 if (b>a)
290 {
291 k1+=m2;
292 k2+=m1;
293 }
294 else
295 {
296 k1+=m1;
297 k2+=m2;
298 }
299 break;
300 }
301 default:
302 {
303 k1=-1;
304 k2=-1;
305 }
306 }
307 }
308 //printf("m1=%d, m2=%d\n",m1,m2);
309 return(k1*N+k2);
310 }
311}
312#endif
313
314static inline int index_sparse_to_full_2d(nsfft_plan *ths, int k)
315{
316 /* only by lookup table */
317 if( k < ths->N_total)
318 return ths->index_sparse_to_full[k];
319 else
320 return -1;
321}
322
323#ifndef NSFTT_DISABLE_TEST
324static int index_full_to_sparse_2d(int J, int k)
325{
326 int N=X(exp2i)(J+2); /* number of full coeffs */
327 int N_B=X(exp2i)(J); /* number in each sparse block */
328
329 int k1=k/N-N/2; /* coordinates in the full grid */
330 int k2=k%N-N/2; /* k1: row, k2: column */
331
332 int r,a,b;
333
334 a=X(exp2i)(J/2+1);
335
336 if ( (k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) )
337 {
338 return(4*((J+1)/2+1)*N_B+(k1+a/2)*a+(k2+a/2));
339 }
340
341 for (r=0; r<=(J+1)/2; r++)
342 {
343 b=X(exp2i)(r);
344 a=X(exp2i)(J-r);
345 if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=a) && (k2<2*a) )
346 {
347 if (a>=b)
348 return((4*r+0)*N_B+(k1+b/2)*a+(k2-a));
349 else
350 return((4*r+0)*N_B+(k2-a)*b+(k1+b/2));
351 }
352 else if ( (k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
353 {
354 if (a>b)
355 return((4*r+1)*N_B+(k2+b/2)*a+(k1-a));
356 else
357 return((4*r+1)*N_B+(k1-a)*b+(k2+b/2));
358 }
359 else if ( (k1>=-(b/2)) && (k1<(b+1)/2) && (k2>=-2*a) && (k2<-a) )
360 {
361 if (a>=b)
362 return((4*r+2)*N_B+(k1+b/2)*a+(k2+2*a));
363 else
364 return((4*r+2)*N_B+(k2+2*a)*b+(k1+b/2));
365 }
366 else if ( (k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) )
367 {
368 if (a>b)
369 return((4*r+3)*N_B+(k2+b/2)*a+(k1+2*a));
370 else
371 return((4*r+3)*N_B+(k1+2*a)*b+(k2+b/2));
372 }
373 }
374
375 return(-1);
376}
377#endif
378
379#ifdef GAUSSIAN
380static void init_index_sparse_to_full_2d(nsfft_plan *ths)
381{
382 int k_S;
383
384 for (k_S=0; k_S<ths->N_total; k_S++)
385 ths->index_sparse_to_full[k_S]=index_sparse_to_full_direct_2d(ths->J, k_S);
386}
387#endif
388
389#ifdef GAUSSIAN
390static inline int index_sparse_to_full_3d(nsfft_plan *ths, int k)
391{
392 /* only by lookup table */
393 if( k < ths->N_total)
394 return ths->index_sparse_to_full[k];
395 else
396 return -1;
397}
398#endif
399
400#ifndef NSFTT_DISABLE_TEST
401static int index_full_to_sparse_3d(int J, int k)
402{
403 int N=X(exp2i)(J+2); /* length of the full grid */
404 int N_B_r; /* size of a sparse block in level r */
405 int sum_N_B_less_r; /* sum N_B_r */
406
407 int r,a,b;
408
409 int k3=(k%N)-N/2; /* coordinates in the full grid */
410 int k2=((k/N)%N)-N/2;
411 int k1=k/(N*N)-N/2;
412
413 a=X(exp2i)(J/2+1); /* length of center block */
414
415 if((k1>=-(a/2)) && (k1<a/2) && (k2>=(-a/2)) && (k2<a/2) && (k3>=(-a/2)) &&
416 (k3<a/2))
417 {
418 return(6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+((k1+a/2)*a+(k2+a/2))*a+
419 (k3+a/2));
420 }
421
422 sum_N_B_less_r=0;
423 for (r=0; r<=(J+1)/2; r++)
424 {
425 a=X(exp2i)(J-r);
426 b=X(exp2i)(r);
427
428 N_B_r=a*b*b;
429
430 /* right - rear - top - left - front - bottom */
431 if ((k1>=a) && (k1<2*a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
432 (k3>=-(b/2)) && (k3<(b+1)/2)) /* right */
433 {
434 if(a>b)
435 return sum_N_B_less_r+N_B_r*0 + ((k2+b/2)*b+k3+b/2)*a + (k1-a);
436 else
437 return sum_N_B_less_r+N_B_r*0 + ((k1-a)*b+(k2+b/2))*b + (k3+b/2);
438 }
439 else if ((k2>=a) && (k2<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
440 (k3>=-(b/2)) && (k3<(b+1)/2)) /* rear */
441 {
442 if(a>b)
443 return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+k3+b/2)*a + (k2-a);
444 else if (a==b)
445 return sum_N_B_less_r+N_B_r*1 + ((k1+b/2)*b+(k2-a))*a + (k3+b/2);
446 else
447 return sum_N_B_less_r+N_B_r*1 + ((k2-a)*b+(k1+b/2))*b + (k3+b/2);
448 }
449 else if ((k3>=a) && (k3<2*a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
450 (k2>=-(b/2)) && (k2<(b+1)/2)) /* top */
451 {
452 if(a>=b)
453 return sum_N_B_less_r+N_B_r*2 + ((k1+b/2)*b+k2+b/2)*a + (k3-a);
454 else
455 return sum_N_B_less_r+N_B_r*2 + ((k3-a)*b+(k1+b/2))*b + (k2+b/2);
456 }
457
458 else if ((k1>=-2*a) && (k1<-a) && (k2>=-(b/2)) && (k2<(b+1)/2) &&
459 (k3>=-(b/2)) && (k3<(b+1)/2)) /* left */
460 {
461 if(a>b)
462 return sum_N_B_less_r+N_B_r*3 + ((k2+b/2)*b+k3+b/2)*a + (k1+2*a);
463 else
464 return sum_N_B_less_r+N_B_r*3 + ((k1+2*a)*b+(k2+b/2))*b + (k3+b/2);
465 }
466 else if ((k2>=-2*a) && (k2<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
467 (k3>=-(b/2)) && (k3<(b+1)/2)) /* front */
468 {
469 if(a>b)
470 return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+k3+b/2)*a + (k2+2*a);
471 else if (a==b)
472 return sum_N_B_less_r+N_B_r*4 + ((k1+b/2)*b+(k2+2*a))*a + (k3+b/2);
473 else
474 return sum_N_B_less_r+N_B_r*4 + ((k2+2*a)*b+(k1+b/2))*b + (k3+b/2);
475 }
476 else if ((k3>=-2*a) && (k3<-a) && (k1>=-(b/2)) && (k1<(b+1)/2) &&
477 (k2>=-(b/2)) && (k2<(b+1)/2)) /* bottom */
478 {
479 if(a>=b)
480 return sum_N_B_less_r+N_B_r*5 + ((k1+b/2)*b+k2+b/2)*a + (k3+2*a);
481 else
482 return sum_N_B_less_r+N_B_r*5 + ((k3+2*a)*b+(k1+b/2))*b + (k2+b/2);
483 }
484
485 sum_N_B_less_r+=6*N_B_r;
486 } /* for(r) */
487
488 return(-1);
489}
490#endif
491
492#ifdef GAUSSIAN
493static void init_index_sparse_to_full_3d(nsfft_plan *ths)
494{
495 int k1,k2,k3,k_s,r;
496 int a,b;
497 int N=X(exp2i)(ths->J+2); /* length of the full grid */
498 int Nc=ths->center_nfft_plan->N[0]; /* length of the center block */
499
500 for (k_s=0, r=0; r<=(ths->J+1)/2; r++)
501 {
502 a=X(exp2i)(ths->J-r);
503 b=X(exp2i)(r);
504
505 /* right - rear - top - left - front - bottom */
506
507 /* right */
508 if(a>b)
509 for(k2=-b/2;k2<(b+1)/2;k2++)
510 for(k3=-b/2;k3<(b+1)/2;k3++)
511 for(k1=a; k1<2*a; k1++,k_s++)
512 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
513 else
514 for(k1=a; k1<2*a; k1++)
515 for(k2=-b/2;k2<(b+1)/2;k2++)
516 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
517 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
518
519 /* rear */
520 if(a>b)
521 for(k1=-b/2;k1<(b+1)/2;k1++)
522 for(k3=-b/2;k3<(b+1)/2;k3++)
523 for(k2=a; k2<2*a; k2++,k_s++)
524 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
525 else if(a==b)
526 for(k1=-b/2;k1<(b+1)/2;k1++)
527 for(k2=a; k2<2*a; k2++)
528 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
529 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
530 else
531 for(k2=a; k2<2*a; k2++)
532 for(k1=-b/2;k1<(b+1)/2;k1++)
533 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
534 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
535
536 /* top */
537 if(a>=b)
538 for(k1=-b/2;k1<(b+1)/2;k1++)
539 for(k2=-b/2;k2<(b+1)/2;k2++)
540 for(k3=a; k3<2*a; k3++,k_s++)
541 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
542 else
543 for(k3=a; k3<2*a; k3++)
544 for(k1=-b/2;k1<(b+1)/2;k1++)
545 for(k2=-b/2;k2<(b+1)/2;k2++,k_s++)
546 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
547
548 /* left */
549 if(a>b)
550 for(k2=-b/2;k2<(b+1)/2;k2++)
551 for(k3=-b/2;k3<(b+1)/2;k3++)
552 for(k1=-2*a; k1<-a; k1++,k_s++)
553 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
554 else
555 for(k1=-2*a; k1<-a; k1++)
556 for(k2=-b/2;k2<(b+1)/2;k2++)
557 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
558 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
559
560 /* front */
561 if(a>b)
562 for(k1=-b/2;k1<(b+1)/2;k1++)
563 for(k3=-b/2;k3<(b+1)/2;k3++)
564 for(k2=-2*a; k2<-a; k2++,k_s++)
565 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
566 else if(a==b)
567 for(k1=-b/2;k1<(b+1)/2;k1++)
568 for(k2=-2*a; k2<-a; k2++)
569 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
570 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
571 else
572 for(k2=-2*a; k2<-a; k2++)
573 for(k1=-b/2;k1<(b+1)/2;k1++)
574 for(k3=-b/2;k3<(b+1)/2;k3++,k_s++)
575 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
576
577 /* top */
578 if(a>=b)
579 for(k1=-b/2;k1<(b+1)/2;k1++)
580 for(k2=-b/2;k2<(b+1)/2;k2++)
581 for(k3=-2*a; k3<-a; k3++,k_s++)
582 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
583 else
584 for(k3=-2*a; k3<-a; k3++)
585 for(k1=-b/2;k1<(b+1)/2;k1++)
586 for(k2=-b/2;k2<(b+1)/2;k2++,k_s++)
587 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
588 }
589
590 /* center */
591 for(k1=-Nc/2;k1<Nc/2;k1++)
592 for(k2=-Nc/2;k2<Nc/2;k2++)
593 for(k3=-Nc/2; k3<Nc/2; k3++,k_s++)
594 ths->index_sparse_to_full[k_s]=((k1+N/2)*N+k2+N/2)*N+k3+N/2;
595}
596#endif
597
598/* copies ths->f_hat to ths_plan->f_hat */
599void nsfft_cp(nsfft_plan *ths, nfft_plan *ths_full_plan)
600{
601 int k;
602
603 /* initialize f_hat with zero values */
604 memset(ths_full_plan->f_hat, 0, ths_full_plan->N_total*sizeof(double _Complex));
605
606 /* copy values at hyperbolic grid points */
607 for(k=0;k<ths->N_total;k++)
608 ths_full_plan->f_hat[ths->index_sparse_to_full[k]]=ths->f_hat[k];
609
610 /* copy nodes */
611 memcpy(ths_full_plan->x,ths->act_nfft_plan->x,ths->M_total*ths->d*sizeof(double));
612}
613
614#ifndef NSFTT_DISABLE_TEST
615/* test copy_sparse_to_full */
616static void test_copy_sparse_to_full_2d(nsfft_plan *ths, nfft_plan *ths_full_plan)
617{
618 int r;
619 int k1, k2;
620 int a,b;
621 const int J=ths->J; /* N=2^J */
622 const int N=ths_full_plan->N[0]; /* size of full NFFT */
623 const int N_B=X(exp2i)(J); /* size of small blocks */
624
625 /* copy sparse plan to full plan */
626 nsfft_cp(ths, ths_full_plan);
627
628 /* show blockwise f_hat */
629 printf("f_hat blockwise\n");
630 for (r=0; r<=(J+1)/2; r++){
631 a=X(exp2i)(J-r); b=X(exp2i)(r);
632
633 printf("top\n");
634 for (k1=0; k1<a; k1++){
635 for (k2=0; k2<b; k2++){
636 printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+1)*N_B+ k1*b+k2]),
637 cimag(ths->f_hat[(4*r+1)*N_B+ k1*b+k2]));
638 }
639 printf("\n");
640 }
641
642 printf("bottom\n");
643 for (k1=0; k1<a; k1++){
644 for (k2=0; k2<b; k2++){
645 printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+3)*N_B+ k1*b+k2]),
646 cimag(ths->f_hat[(4*r+3)*N_B+ k1*b+k2]));
647 }
648 printf("\n");
649 }
650
651 printf("right\n");
652 for (k2=0; k2<b; k2++){
653 for (k1=0; k1<a; k1++){
654 printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+0)*N_B+ k2*a+k1]),
655 cimag(ths->f_hat[(4*r+0)*N_B+ k2*a+k1]));
656 }
657 printf("\n");
658 }
659
660 printf("left\n");
661 for (k2=0; k2<b; k2++){
662 for (k1=0; k1<a; k1++){
663 printf("(%1.1f,%1.1f) ", creal(ths->f_hat[(4*r+2)*N_B+ k2*a+k1]),
664 cimag(ths->f_hat[(4*r+2)*N_B+ k2*a+k1]));
665 }
666 printf("\n");
667 }
668 }
669
670 return;
671 /* show full f_hat */
672 printf("full f_hat\n");
673 for (k1=0;k1<N;k1++){
674 for (k2=0;k2<N;k2++){
675 printf("(%1.1f,%1.1f) ", creal(ths_full_plan->f_hat[k1*N+k2]),
676 cimag(ths_full_plan->f_hat[k1*N+k2]));
677 }
678 printf("\n");
679 }
680}
681#endif
682
683#ifndef NSFTT_DISABLE_TEST
684static void test_sparse_to_full_2d(nsfft_plan* ths)
685{
686 int k_S,k1,k2;
687 int N=X(exp2i)(ths->J+2);
688
689 printf("N=%d\n\n",N);
690
691 for(k1=0;k1<N;k1++)
692 for(k2=0;k2<N;k2++)
693 {
694 k_S=index_full_to_sparse_2d(ths->J, k1*N+k2);
695 if(k_S!=-1)
696 printf("(%+d, %+d)\t= %+d \t= %+d = %+d \n",k1-N/2,k2-N/2,
697 k1*N+k2, k_S, ths->index_sparse_to_full[k_S]);
698 }
699}
700#endif
701
702#ifndef NSFTT_DISABLE_TEST
703static void test_sparse_to_full_3d(nsfft_plan* ths)
704{
705 int k_S,k1,k2,k3;
706 int N=X(exp2i)(ths->J+2);
707
708 printf("N=%d\n\n",N);
709
710 for(k1=0;k1<N;k1++)
711 for(k2=0;k2<N;k2++)
712 for(k3=0;k3<N;k3++)
713 {
714 k_S=index_full_to_sparse_3d(ths->J, (k1*N+k2)*N+k3);
715 if(k_S!=-1)
716 printf("(%d, %d, %d)\t= %d \t= %d = %d \n",k1-N/2,k2-N/2,k3-N/2,
717 (k1*N+k2)*N+k3,k_S, ths->index_sparse_to_full[k_S]);
718 }
719}
720#endif
721
722
723void nsfft_init_random_nodes_coeffs(nsfft_plan *ths)
724{
725 int j;
726
727 /* init frequencies */
729
730 /* init nodes */
732
733 if(ths->d==2)
734 for(j=0;j<ths->M_total;j++)
735 {
736 ths->x_transposed[2*j+0]=ths->act_nfft_plan->x[2*j+1];
737 ths->x_transposed[2*j+1]=ths->act_nfft_plan->x[2*j+0];
738 }
739 else /* this->d==3 */
740 for(j=0;j<ths->M_total;j++)
741 {
742 ths->x_102[3*j+0]=ths->act_nfft_plan->x[3*j+1];
743 ths->x_102[3*j+1]=ths->act_nfft_plan->x[3*j+0];
744 ths->x_102[3*j+2]=ths->act_nfft_plan->x[3*j+2];
745
746 ths->x_201[3*j+0]=ths->act_nfft_plan->x[3*j+2];
747 ths->x_201[3*j+1]=ths->act_nfft_plan->x[3*j+0];
748 ths->x_201[3*j+2]=ths->act_nfft_plan->x[3*j+1];
749
750 ths->x_120[3*j+0]=ths->act_nfft_plan->x[3*j+1];
751 ths->x_120[3*j+1]=ths->act_nfft_plan->x[3*j+2];
752 ths->x_120[3*j+2]=ths->act_nfft_plan->x[3*j+0];
753
754 ths->x_021[3*j+0]=ths->act_nfft_plan->x[3*j+0];
755 ths->x_021[3*j+1]=ths->act_nfft_plan->x[3*j+2];
756 ths->x_021[3*j+2]=ths->act_nfft_plan->x[3*j+1];
757 }
758}
759
760static void nsdft_trafo_2d(nsfft_plan *ths)
761{
762 int j,k_S,k_L,k0,k1;
763 double omega;
764 int N=X(exp2i)(ths->J+2);
765
766 memset(ths->f,0,ths->M_total*sizeof(double _Complex));
767
768 for(k_S=0;k_S<ths->N_total;k_S++)
769 {
770 k_L=ths->index_sparse_to_full[k_S];
771 k0=k_L / N;
772 k1=k_L % N;
773
774 for(j=0;j<ths->M_total;j++)
775 {
776 omega =
777 ((double)(k0 - N/2)) * ths->act_nfft_plan->x[2 * j + 0] +
778 ((double)(k1 - N/2)) * ths->act_nfft_plan->x[2 * j + 1];
779 ths->f[j] += ths->f_hat[k_S] * cexp( - I*2*KPI*omega);
780 }
781 }
782} /* void nsdft_trafo_2d */
783
784static void nsdft_trafo_3d(nsfft_plan *ths)
785{
786 int j,k_S,k0,k1,k2;
787 double omega;
788 int N=X(exp2i)(ths->J+2);
789 int k_L;
790
791 memset(ths->f,0,ths->M_total*sizeof(double _Complex));
792
793 for(k_S=0;k_S<ths->N_total;k_S++)
794 {
795 k_L=ths->index_sparse_to_full[k_S];
796
797 k0=k_L/(N*N);
798 k1=(k_L/N)%N;
799 k2=k_L%N;
800
801 for(j=0;j<ths->M_total;j++)
802 {
803 omega =
804 ((double)(k0 - N/2)) * ths->act_nfft_plan->x[3 * j + 0] +
805 ((double)(k1 - N/2)) * ths->act_nfft_plan->x[3 * j + 1] +
806 ((double)(k2 - N/2)) * ths->act_nfft_plan->x[3 * j + 2];
807 ths->f[j] += ths->f_hat[k_S] * cexp( - I*2*KPI*omega);
808 }
809 }
810} /* void nsdft_trafo_3d */
811
812void nsfft_trafo_direct(nsfft_plan *ths)
813{
814 if(ths->d==2)
815 nsdft_trafo_2d(ths);
816 else
817 nsdft_trafo_3d(ths);
818}
819
820static void nsdft_adjoint_2d(nsfft_plan *ths)
821{
822 int j,k_S,k_L,k0,k1;
823 double omega;
824 int N=X(exp2i)(ths->J+2);
825
826 memset(ths->f_hat,0,ths->N_total*sizeof(double _Complex));
827
828 for(k_S=0;k_S<ths->N_total;k_S++)
829 {
830 k_L=ths->index_sparse_to_full[k_S];
831 k0=k_L / N;
832 k1=k_L % N;
833
834 for(j=0;j<ths->M_total;j++)
835 {
836 omega =
837 ((double)(k0 - N/2)) * ths->act_nfft_plan->x[2 * j + 0] +
838 ((double)(k1 - N/2)) * ths->act_nfft_plan->x[2 * j + 1];
839 ths->f_hat[k_S] += ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
840 }
841 }
842} /* void nsdft_adjoint_2d */
843
844static void nsdft_adjoint_3d(nsfft_plan *ths)
845{
846 int j,k_S,k0,k1,k2;
847 double omega;
848 int N=X(exp2i)(ths->J+2);
849 int k_L;
850
851 memset(ths->f_hat,0,ths->N_total*sizeof(double _Complex));
852
853 for(k_S=0;k_S<ths->N_total;k_S++)
854 {
855 k_L=ths->index_sparse_to_full[k_S];
856
857 k0=k_L/(N*N);
858 k1=(k_L/N)%N;
859 k2=k_L%N;
860
861 for(j=0;j<ths->M_total;j++)
862 {
863 omega =
864 ((double)(k0 - N/2)) * ths->act_nfft_plan->x[3 * j + 0] +
865 ((double)(k1 - N/2)) * ths->act_nfft_plan->x[3 * j + 1] +
866 ((double)(k2 - N/2)) * ths->act_nfft_plan->x[3 * j + 2];
867 ths->f_hat[k_S] += ths->f[j] * cexp( + _Complex_I*2*KPI*omega);
868 }
869 }
870} /* void nsdft_adjoint_3d */
871
872void nsfft_adjoint_direct(nsfft_plan *ths)
873{
874 if(ths->d==2)
875 nsdft_adjoint_2d(ths);
876 else
877 nsdft_adjoint_3d(ths);
878}
879
880static void nsfft_trafo_2d(nsfft_plan *ths)
881{
882 int r,rr,j;
883 double temp;
884
885 int M=ths->M_total;
886 int J=ths->J;
887
888 /* center */
889 ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*X(exp2i)(J);
890
891 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
892 nfft_trafo_direct(ths->center_nfft_plan);
893 else
894 nfft_trafo(ths->center_nfft_plan);
895
896 for (j=0; j<M; j++)
897 ths->f[j] = ths->center_nfft_plan->f[j];
898
899 for(rr=0;rr<=(J+1)/2;rr++)
900 {
901 r=MIN(rr,J-rr);
902 ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[r];
903 ths->act_nfft_plan->N[0]=X(exp2i)(r); ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
904 ths->act_nfft_plan->N[1]=X(exp2i)(J-r); ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
905
906 /*printf("%d x %d\n",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1]);*/
907
908 temp=-3.0*KPI*X(exp2i)(J-rr);
909
910 /* right */
911 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*X(exp2i)(J);
912
913 if(r<rr)
914 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
915
916 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
917 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
918 nfft_trafo_direct(ths->act_nfft_plan);
919 else
920 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
921 else
922 nfft_trafo(ths->act_nfft_plan);
923
924 if(r<rr)
925 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
926
927 for (j=0; j<M; j++)
928 ths->f[j] += ths->act_nfft_plan->f[j] *
929 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+1]);
930
931 /* top */
932 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*X(exp2i)(J);
933
934 if((r==rr)&&(J-rr!=rr))
935 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
936
937 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
938 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
939 nfft_trafo_direct(ths->act_nfft_plan);
940 else
941 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
942 else
943 nfft_trafo(ths->act_nfft_plan);
944
945 if((r==rr)&&(J-rr!=rr))
946 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
947
948 for (j=0; j<M; j++)
949 ths->f[j] += ths->act_nfft_plan->f[j] *
950 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+0]);
951
952 /* left */
953 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*X(exp2i)(J);
954
955 if(r<rr)
956 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
957
958 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
959 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
960 nfft_trafo_direct(ths->act_nfft_plan);
961 else
962 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
963 else
964 nfft_trafo(ths->act_nfft_plan);
965
966 if(r<rr)
967 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
968
969 for (j=0; j<M; j++)
970 ths->f[j] += ths->act_nfft_plan->f[j] *
971 cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
972
973 /* bottom */
974 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*X(exp2i)(J);
975
976 if((r==rr)&&(J-rr!=rr))
977 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
978
979 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
980 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
981 nfft_trafo_direct(ths->act_nfft_plan);
982 else
983 short_nfft_trafo_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
984 else
985 nfft_trafo(ths->act_nfft_plan);
986
987 if((r==rr)&&(J-rr!=rr))
988 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
989
990 for (j=0; j<M; j++)
991 ths->f[j] += ths->act_nfft_plan->f[j] *
992 cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
993 } /* for(rr) */
994} /* void nsfft_trafo_2d */
995
996static void nsfft_adjoint_2d(nsfft_plan *ths)
997{
998 int r,rr,j;
999 double temp;
1000
1001 int M=ths->M_total;
1002 int J=ths->J;
1003
1004 /* center */
1005 for (j=0; j<M; j++)
1006 ths->center_nfft_plan->f[j] = ths->f[j];
1007
1008 ths->center_nfft_plan->f_hat=ths->f_hat+4*((J+1)/2+1)*X(exp2i)(J);
1009
1010 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1011 nfft_adjoint_direct(ths->center_nfft_plan);
1012 else
1013 nfft_adjoint(ths->center_nfft_plan);
1014
1015 for(rr=0;rr<=(J+1)/2;rr++)
1016 {
1017 r=MIN(rr,J-rr);
1018 ths->act_nfft_plan->my_fftw_plan2 = ths->set_fftw_plan2[r];
1019 ths->act_nfft_plan->N[0]=X(exp2i)(r); ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
1020 ths->act_nfft_plan->N[1]=X(exp2i)(J-r); ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
1021
1022 /*printf("%d x %d\n",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1]);*/
1023
1024 temp=-3.0*KPI*X(exp2i)(J-rr);
1025
1026 /* right */
1027 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+0)*X(exp2i)(J);
1028
1029 for (j=0; j<M; j++)
1030 ths->act_nfft_plan->f[j]= ths->f[j] *
1031 cexp( - I*temp*ths->act_nfft_plan->x[2*j+1]);
1032
1033 if(r<rr)
1034 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1035
1036 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1037 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1038 nfft_adjoint_direct(ths->act_nfft_plan);
1039 else
1040 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1041 else
1042 nfft_adjoint(ths->act_nfft_plan);
1043
1044 if(r<rr)
1045 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1046
1047 /* top */
1048 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+1)*X(exp2i)(J);
1049
1050 for (j=0; j<M; j++)
1051 ths->act_nfft_plan->f[j]= ths->f[j] *
1052 cexp( - I*temp*ths->act_nfft_plan->x[2*j+0]);
1053
1054 if((r==rr)&&(J-rr!=rr))
1055 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1056
1057 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1058 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1059 nfft_adjoint_direct(ths->act_nfft_plan);
1060 else
1061 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1062 else
1063 nfft_adjoint(ths->act_nfft_plan);
1064
1065 if((r==rr)&&(J-rr!=rr))
1066 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1067
1068 /* left */
1069 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+2)*X(exp2i)(J);
1070
1071 for (j=0; j<M; j++)
1072 ths->act_nfft_plan->f[j]= ths->f[j] *
1073 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+1]);
1074
1075 if(r<rr)
1076 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1077
1078 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1079 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1080 nfft_adjoint_direct(ths->act_nfft_plan);
1081 else
1082 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1083 else
1084 nfft_adjoint(ths->act_nfft_plan);
1085
1086 if(r<rr)
1087 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1088
1089 /* bottom */
1090 ths->act_nfft_plan->f_hat=ths->f_hat+(4*rr+3)*X(exp2i)(J);
1091
1092 for (j=0; j<M; j++)
1093 ths->act_nfft_plan->f[j]= ths->f[j] *
1094 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[2*j+0]);
1095
1096 if((r==rr)&&(J-rr!=rr))
1097 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1098
1099 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1100 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1101 nfft_adjoint_direct(ths->act_nfft_plan);
1102 else
1103 short_nfft_adjoint_2d(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1104 else
1105 nfft_adjoint(ths->act_nfft_plan);
1106
1107 if((r==rr)&&(J-rr!=rr))
1108 RSWAP(ths->act_nfft_plan->x,ths->x_transposed);
1109 } /* for(rr) */
1110} /* void nsfft_adjoint_2d */
1111
1112static void nsfft_trafo_3d(nsfft_plan *ths)
1113{
1114 int r,rr,j;
1115 double temp;
1116 int sum_N_B_less_r,N_B_r,a,b;
1117
1118 int M=ths->M_total;
1119 int J=ths->J;
1120
1121 /* center */
1122 ths->center_nfft_plan->f_hat=ths->f_hat+6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1);
1123
1124 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1125 nfft_trafo_direct(ths->center_nfft_plan);
1126 else
1127 nfft_trafo(ths->center_nfft_plan);
1128
1129 for (j=0; j<M; j++)
1130 ths->f[j] = ths->center_nfft_plan->f[j];
1131
1132 sum_N_B_less_r=0;
1133 for(rr=0;rr<=(J+1)/2;rr++)
1134 {
1135 a=X(exp2i)(J-rr);
1136 b=X(exp2i)(rr);
1137
1138 N_B_r=a*b*b;
1139
1140 r=MIN(rr,J-rr);
1141 ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
1142
1143 ths->act_nfft_plan->N[0]=X(exp2i)(r);
1144 if(a<b)
1145 ths->act_nfft_plan->N[1]=X(exp2i)(J-r);
1146 else
1147 ths->act_nfft_plan->N[1]=X(exp2i)(r);
1148 ths->act_nfft_plan->N[2]=X(exp2i)(J-r);
1149
1150 /*printf("\n\n%d x %d x %d:\t",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1],ths->act_nfft_plan->N[2]); fflush(stdout);*/
1151
1152 ths->act_nfft_plan->N_total=ths->act_nfft_plan->N[0]*ths->act_nfft_plan->N[1]*ths->act_nfft_plan->N[2];
1153 ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
1154 ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
1155 ths->act_nfft_plan->n[2]=ths->sigma*ths->act_nfft_plan->N[2];
1156 ths->act_nfft_plan->n_total=ths->act_nfft_plan->n[0]*ths->act_nfft_plan->n[1]*ths->act_nfft_plan->n[2];
1157
1158 /* only for right - rear - top */
1159 if((J==0)||((J==1)&&(rr==1)))
1160 temp=-2.0*KPI;
1161 else
1162 temp=-3.0*KPI*X(exp2i)(J-rr);
1163
1164 /* right */
1165 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
1166
1167 if(a>b)
1168 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1169
1170 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1171 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1172 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1173 nfft_trafo_direct(ths->act_nfft_plan);
1174 else
1175 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1176 else
1177 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1178 else
1179 nfft_trafo(ths->act_nfft_plan);
1180
1181 if(a>b)
1182 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1183
1184 for (j=0; j<M; j++)
1185 ths->f[j] += ths->act_nfft_plan->f[j] *
1186 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+0]);
1187
1188 /* rear */
1189 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
1190
1191 if(a>b)
1192 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1193 if(a<b)
1194 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1195
1196 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1197 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1198 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1199 nfft_trafo_direct(ths->act_nfft_plan);
1200 else
1201 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1202 else
1203 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1204 else
1205 nfft_trafo(ths->act_nfft_plan);
1206
1207 if(a>b)
1208 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1209 if(a<b)
1210 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1211
1212 for (j=0; j<M; j++)
1213 ths->f[j] += ths->act_nfft_plan->f[j] *
1214 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+1]);
1215
1216 /* top */
1217 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
1218
1219 if(a<b)
1220 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1221
1222 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1223 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1224 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1225 nfft_trafo_direct(ths->act_nfft_plan);
1226 else
1227 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1228 else
1229 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1230 else
1231 nfft_trafo(ths->act_nfft_plan);
1232
1233 if(a<b)
1234 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1235
1236 for (j=0; j<M; j++)
1237 ths->f[j] += ths->act_nfft_plan->f[j] *
1238 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+2]);
1239
1240 /* only for left - front - bottom */
1241 if((J==0)||((J==1)&&(rr==1)))
1242 temp=-4.0*KPI;
1243 else
1244 temp=-3.0*KPI*X(exp2i)(J-rr);
1245
1246 /* left */
1247 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
1248
1249 if(a>b)
1250 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1251
1252 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1253 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1254 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1255 nfft_trafo_direct(ths->act_nfft_plan);
1256 else
1257 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1258 else
1259 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1260 else
1261 nfft_trafo(ths->act_nfft_plan);
1262
1263 if(a>b)
1264 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1265
1266 for (j=0; j<M; j++)
1267 ths->f[j] += ths->act_nfft_plan->f[j] *
1268 cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
1269
1270 /* front */
1271 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
1272
1273 if(a>b)
1274 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1275 if(a<b)
1276 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1277
1278 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1279 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1280 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1281 nfft_trafo_direct(ths->act_nfft_plan);
1282 else
1283 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1284 else
1285 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1286 else
1287 nfft_trafo(ths->act_nfft_plan);
1288
1289 if(a>b)
1290 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1291 if(a<b)
1292 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1293
1294 for (j=0; j<M; j++)
1295 ths->f[j] += ths->act_nfft_plan->f[j] *
1296 cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
1297
1298 /* bottom */
1299 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
1300
1301 if(a<b)
1302 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1303
1304 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1305 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1306 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1307 nfft_trafo_direct(ths->act_nfft_plan);
1308 else
1309 short_nfft_trafo_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1310 else
1311 short_nfft_trafo_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1312 else
1313 nfft_trafo(ths->act_nfft_plan);
1314
1315 if(a<b)
1316 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1317
1318 for (j=0; j<M; j++)
1319 ths->f[j] += ths->act_nfft_plan->f[j] *
1320 cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
1321
1322 sum_N_B_less_r+=6*N_B_r;
1323 } /* for(rr) */
1324} /* void nsfft_trafo_3d */
1325
1326static void nsfft_adjoint_3d(nsfft_plan *ths)
1327{
1328 int r,rr,j;
1329 double temp;
1330 int sum_N_B_less_r,N_B_r,a,b;
1331
1332 int M=ths->M_total;
1333 int J=ths->J;
1334
1335 /* center */
1336 for (j=0; j<M; j++)
1337 ths->center_nfft_plan->f[j] = ths->f[j];
1338
1339 ths->center_nfft_plan->f_hat=ths->f_hat+6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1);
1340
1341 if (ths->center_nfft_plan->N[0]<=ths->center_nfft_plan->m)
1342 nfft_adjoint_direct(ths->center_nfft_plan);
1343 else
1344 nfft_adjoint(ths->center_nfft_plan);
1345
1346 sum_N_B_less_r=0;
1347 for(rr=0;rr<=(J+1)/2;rr++)
1348 {
1349 a=X(exp2i)(J-rr);
1350 b=X(exp2i)(rr);
1351
1352 N_B_r=a*b*b;
1353
1354 r=MIN(rr,J-rr);
1355 ths->act_nfft_plan->my_fftw_plan1 = ths->set_fftw_plan1[rr];
1356 ths->act_nfft_plan->my_fftw_plan2 = ths->set_fftw_plan2[rr];
1357
1358 ths->act_nfft_plan->N[0]=X(exp2i)(r);
1359 if(a<b)
1360 ths->act_nfft_plan->N[1]=X(exp2i)(J-r);
1361 else
1362 ths->act_nfft_plan->N[1]=X(exp2i)(r);
1363 ths->act_nfft_plan->N[2]=X(exp2i)(J-r);
1364
1365 /*printf("\n\n%d x %d x %d:\t",ths->act_nfft_plan->N[0],ths->act_nfft_plan->N[1],ths->act_nfft_plan->N[2]); fflush(stdout);*/
1366
1367 ths->act_nfft_plan->N_total=ths->act_nfft_plan->N[0]*ths->act_nfft_plan->N[1]*ths->act_nfft_plan->N[2];
1368 ths->act_nfft_plan->n[0]=ths->sigma*ths->act_nfft_plan->N[0];
1369 ths->act_nfft_plan->n[1]=ths->sigma*ths->act_nfft_plan->N[1];
1370 ths->act_nfft_plan->n[2]=ths->sigma*ths->act_nfft_plan->N[2];
1371 ths->act_nfft_plan->n_total=ths->act_nfft_plan->n[0]*ths->act_nfft_plan->n[1]*ths->act_nfft_plan->n[2];
1372
1373 /* only for right - rear - top */
1374 if((J==0)||((J==1)&&(rr==1)))
1375 temp=-2.0*KPI;
1376 else
1377 temp=-3.0*KPI*X(exp2i)(J-rr);
1378
1379 /* right */
1380 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*0;
1381
1382 for (j=0; j<M; j++)
1383 ths->act_nfft_plan->f[j]= ths->f[j] *
1384 cexp( - I*temp*ths->act_nfft_plan->x[3*j+0]);
1385
1386 if(a>b)
1387 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1388
1389 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1390 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1391 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1392 nfft_adjoint_direct(ths->act_nfft_plan);
1393 else
1394 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1395 else
1396 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1397 else
1398 nfft_adjoint(ths->act_nfft_plan);
1399
1400 if(a>b)
1401 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1402
1403 /* rear */
1404 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*1;
1405
1406 for (j=0; j<M; j++)
1407 ths->act_nfft_plan->f[j]= ths->f[j] *
1408 cexp( - I*temp*ths->act_nfft_plan->x[3*j+1]);
1409
1410 if(a>b)
1411 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1412 if(a<b)
1413 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1414
1415 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1416 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1417 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1418 nfft_adjoint_direct(ths->act_nfft_plan);
1419 else
1420 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1421 else
1422 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1423 else
1424 nfft_adjoint(ths->act_nfft_plan);
1425
1426 if(a>b)
1427 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1428 if(a<b)
1429 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1430
1431 /* top */
1432 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*2;
1433
1434 for (j=0; j<M; j++)
1435 ths->act_nfft_plan->f[j]= ths->f[j] *
1436 cexp( - I*temp*ths->act_nfft_plan->x[3*j+2]);
1437
1438 if(a<b)
1439 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1440
1441 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1442 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1443 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1444 nfft_adjoint_direct(ths->act_nfft_plan);
1445 else
1446 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1447 else
1448 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1449 else
1450 nfft_adjoint(ths->act_nfft_plan);
1451
1452 if(a<b)
1453 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1454
1455 /* only for left - front - bottom */
1456 if((J==0)||((J==1)&&(rr==1)))
1457 temp=-4.0*KPI;
1458 else
1459 temp=-3.0*KPI*X(exp2i)(J-rr);
1460
1461 /* left */
1462 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*3;
1463
1464 for (j=0; j<M; j++)
1465 ths->act_nfft_plan->f[j]= ths->f[j] *
1466 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+0]);
1467
1468 if(a>b)
1469 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1470
1471 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1472 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1473 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1474 nfft_adjoint_direct(ths->act_nfft_plan);
1475 else
1476 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1477 else
1478 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1479 else
1480 nfft_adjoint(ths->act_nfft_plan);
1481
1482 if(a>b)
1483 RSWAP(ths->act_nfft_plan->x,ths->x_120);
1484
1485 /* front */
1486 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*4;
1487
1488 for (j=0; j<M; j++)
1489 ths->act_nfft_plan->f[j]= ths->f[j] *
1490 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+1]);
1491
1492 if(a>b)
1493 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1494 if(a<b)
1495 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1496
1497 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1498 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1499 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1500 nfft_adjoint_direct(ths->act_nfft_plan);
1501 else
1502 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1503 else
1504 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1505 else
1506 nfft_adjoint(ths->act_nfft_plan);
1507
1508 if(a>b)
1509 RSWAP(ths->act_nfft_plan->x,ths->x_021);
1510 if(a<b)
1511 RSWAP(ths->act_nfft_plan->x,ths->x_102);
1512
1513 /* bottom */
1514 ths->act_nfft_plan->f_hat=ths->f_hat + sum_N_B_less_r + N_B_r*5;
1515
1516 for (j=0; j<M; j++)
1517 ths->act_nfft_plan->f[j]= ths->f[j] *
1518 cexp( + _Complex_I*temp*ths->act_nfft_plan->x[3*j+2]);
1519
1520 if(a<b)
1521 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1522
1523 if(ths->act_nfft_plan->N[0]<=ths->act_nfft_plan->m)
1524 if(ths->act_nfft_plan->N[1]<=ths->act_nfft_plan->m)
1525 if(ths->act_nfft_plan->N[2]<=ths->act_nfft_plan->m)
1526 nfft_adjoint_direct(ths->act_nfft_plan);
1527 else
1528 short_nfft_adjoint_3d_1(ths->act_nfft_plan,&(ths->set_nfft_plan_1d[r]));
1529 else
1530 short_nfft_adjoint_3d_2(ths->act_nfft_plan,&(ths->set_nfft_plan_2d[r]));
1531 else
1532 nfft_adjoint(ths->act_nfft_plan);
1533
1534 if(a<b)
1535 RSWAP(ths->act_nfft_plan->x,ths->x_201);
1536
1537 sum_N_B_less_r+=6*N_B_r;
1538 } /* for(rr) */
1539} /* void nsfft_adjoint_3d */
1540
1541void nsfft_trafo(nsfft_plan *ths)
1542{
1543 if(ths->d==2)
1544 nsfft_trafo_2d(ths);
1545 else
1546 nsfft_trafo_3d(ths);
1547}
1548
1549void nsfft_adjoint(nsfft_plan *ths)
1550{
1551 if(ths->d==2)
1552 nsfft_adjoint_2d(ths);
1553 else
1554 nsfft_adjoint_3d(ths);
1555}
1556
1557
1558/*========================================================*/
1559/* J >1, no precomputation at all!! */
1560#ifdef GAUSSIAN
1561static void nsfft_init_2d(nsfft_plan *ths, int J, int M, int m, unsigned snfft_flags)
1562{
1563 int r;
1564 int N[2];
1565 int n[2];
1566
1567 ths->flags=snfft_flags;
1568 ths->sigma=2;
1569 ths->J=J;
1570 ths->M_total=M;
1571 ths->N_total=(J+4)*X(exp2i)(J+1);
1572
1573 /* memory allocation */
1574 ths->f = (double _Complex *)nfft_malloc(M*sizeof(double _Complex));
1575 ths->f_hat = (double _Complex *)nfft_malloc(ths->N_total*sizeof(double _Complex));
1576 ths->x_transposed= (double*)nfft_malloc(2*M*sizeof(double));
1577
1578 ths->act_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
1580
1581 ths->set_fftw_plan1=(fftw_plan*) nfft_malloc((J/2+1)*sizeof(fftw_plan));
1582 ths->set_fftw_plan2=(fftw_plan*) nfft_malloc((J/2+1)*sizeof(fftw_plan));
1583
1584 ths->set_nfft_plan_1d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
1585
1586 /* planning the small nffts */
1587 /* r=0 */
1588 N[0]=1; n[0]=ths->sigma*N[0];
1589 N[1]=X(exp2i)(J); n[1]=ths->sigma*N[1];
1590
1591 nfft_init_guru(ths->act_nfft_plan,2,N,M,n,m, FG_PSI| MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
1592
1593 if(ths->act_nfft_plan->flags & PRE_ONE_PSI)
1594 nfft_precompute_one_psi(ths->act_nfft_plan);
1595
1596 ths->set_fftw_plan1[0]=ths->act_nfft_plan->my_fftw_plan1;
1597 ths->set_fftw_plan2[0]=ths->act_nfft_plan->my_fftw_plan2;
1598
1599 for(r=1;r<=J/2;r++)
1600 {
1601 N[0]=X(exp2i)(r); n[0]=ths->sigma*N[0];
1602 N[1]=X(exp2i)(J-r); n[1]=ths->sigma*N[1];
1603 ths->set_fftw_plan1[r] =
1604 fftw_plan_dft(2, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
1605 FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
1606
1607 ths->set_fftw_plan2[r] =
1608 fftw_plan_dft(2, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
1609 FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
1610 }
1611
1612 /* planning the 1d nffts */
1613 for(r=0;r<=X(log2i)(m);r++)
1614 {
1615 N[0]=X(exp2i)(J-r); n[0]=ths->sigma*N[0]; /* ==N[1] of the 2 dimensional plan */
1616
1617 nfft_init_guru(&(ths->set_nfft_plan_1d[r]),1,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
1618 ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags | FG_PSI;
1619 ths->set_nfft_plan_1d[r].K=ths->act_nfft_plan->K;
1620 ths->set_nfft_plan_1d[r].psi=ths->act_nfft_plan->psi;
1621 }
1622
1623 /* center plan */
1624 /* J/2 == floor(((double)J) / 2.0) */
1625 N[0]=X(exp2i)(J/2+1); n[0]=ths->sigma*N[0];
1626 N[1]=X(exp2i)(J/2+1); n[1]=ths->sigma*N[1];
1627 nfft_init_guru(ths->center_nfft_plan,2,N,M,n, m, MALLOC_F| FFTW_INIT,
1628 FFTW_MEASURE);
1629 ths->center_nfft_plan->x= ths->act_nfft_plan->x;
1630 ths->center_nfft_plan->flags = ths->center_nfft_plan->flags|
1631 FG_PSI;
1632 ths->center_nfft_plan->K=ths->act_nfft_plan->K;
1633 ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
1634
1635 if(ths->flags & NSDFT)
1636 {
1637 ths->index_sparse_to_full=(int*)nfft_malloc(ths->N_total*sizeof(int));
1638 init_index_sparse_to_full_2d(ths);
1639 }
1640}
1641#endif
1642
1643/*========================================================*/
1644/* J >1, no precomputation at all!! */
1645#ifdef GAUSSIAN
1646static void nsfft_init_3d(nsfft_plan *ths, int J, int M, int m, unsigned snfft_flags)
1647{
1648 int r,rr,a,b;
1649 int N[3];
1650 int n[3];
1651
1652 ths->flags=snfft_flags;
1653 ths->sigma=2;
1654 ths->J=J;
1655 ths->M_total=M;
1656 ths->N_total=6*X(exp2i)(J)*(X(exp2i)((J+1)/2+1)-1)+X(exp2i)(3*(J/2+1));
1657
1658 /* memory allocation */
1659 ths->f = (double _Complex *)nfft_malloc(M*sizeof(double _Complex));
1660 ths->f_hat = (double _Complex *)nfft_malloc(ths->N_total*sizeof(double _Complex));
1661
1662 ths->x_102= (double*)nfft_malloc(3*M*sizeof(double));
1663 ths->x_201= (double*)nfft_malloc(3*M*sizeof(double));
1664 ths->x_120= (double*)nfft_malloc(3*M*sizeof(double));
1665 ths->x_021= (double*)nfft_malloc(3*M*sizeof(double));
1666
1667 ths->act_nfft_plan = (nfft_plan*)nfft_malloc(sizeof(nfft_plan));
1669
1670 ths->set_fftw_plan1=(fftw_plan*) nfft_malloc(((J+1)/2+1)*sizeof(fftw_plan));
1671 ths->set_fftw_plan2=(fftw_plan*) nfft_malloc(((J+1)/2+1)*sizeof(fftw_plan));
1672
1673 ths->set_nfft_plan_1d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
1674 ths->set_nfft_plan_2d = (nfft_plan*) nfft_malloc((X(log2i)(m)+1)*(sizeof(nfft_plan)));
1675
1676 /* planning the small nffts */
1677 /* r=0 */
1678 N[0]=1; n[0]=ths->sigma*N[0];
1679 N[1]=1; n[1]=ths->sigma*N[1];
1680 N[2]=X(exp2i)(J); n[2]=ths->sigma*N[2];
1681
1682 nfft_init_guru(ths->act_nfft_plan,3,N,M,n,m, FG_PSI| MALLOC_X| MALLOC_F, FFTW_MEASURE);
1683
1684 if(ths->act_nfft_plan->flags & PRE_ONE_PSI)
1685 nfft_precompute_one_psi(ths->act_nfft_plan);
1686
1687 /* malloc g1, g2 for maximal size */
1688 ths->act_nfft_plan->g1 = nfft_malloc(ths->sigma*ths->sigma*ths->sigma*X(exp2i)(J+(J+1)/2)*sizeof(double _Complex));
1689 ths->act_nfft_plan->g2 = nfft_malloc(ths->sigma*ths->sigma*ths->sigma*X(exp2i)(J+(J+1)/2)*sizeof(double _Complex));
1690
1691 ths->act_nfft_plan->my_fftw_plan1 =
1692 fftw_plan_dft(3, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
1693 FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
1694 ths->act_nfft_plan->my_fftw_plan2 =
1695 fftw_plan_dft(3, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
1696 FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
1697
1698 ths->set_fftw_plan1[0]=ths->act_nfft_plan->my_fftw_plan1;
1699 ths->set_fftw_plan2[0]=ths->act_nfft_plan->my_fftw_plan2;
1700
1701 for(rr=1;rr<=(J+1)/2;rr++)
1702 {
1703 a=X(exp2i)(J-rr);
1704 b=X(exp2i)(rr);
1705
1706 r=MIN(rr,J-rr);
1707
1708 n[0]=ths->sigma*X(exp2i)(r);
1709 if(a<b)
1710 n[1]=ths->sigma*X(exp2i)(J-r);
1711 else
1712 n[1]=ths->sigma*X(exp2i)(r);
1713 n[2]=ths->sigma*X(exp2i)(J-r);
1714
1715 ths->set_fftw_plan1[rr] =
1716 fftw_plan_dft(3, n, ths->act_nfft_plan->g1, ths->act_nfft_plan->g2,
1717 FFTW_FORWARD, ths->act_nfft_plan->fftw_flags);
1718 ths->set_fftw_plan2[rr] =
1719 fftw_plan_dft(3, n, ths->act_nfft_plan->g2, ths->act_nfft_plan->g1,
1720 FFTW_BACKWARD, ths->act_nfft_plan->fftw_flags);
1721 }
1722
1723 /* planning the 1d nffts */
1724 for(r=0;r<=X(log2i)(m);r++)
1725 {
1726 N[0]=X(exp2i)(J-r); n[0]=ths->sigma*N[0];
1727 N[1]=X(exp2i)(J-r); n[1]=ths->sigma*N[1];
1728
1729 if(N[0]>m)
1730 {
1731 nfft_init_guru(&(ths->set_nfft_plan_1d[r]),1,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
1732 ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags | FG_PSI;
1733 ths->set_nfft_plan_1d[r].K=ths->act_nfft_plan->K;
1734 ths->set_nfft_plan_1d[r].psi=ths->act_nfft_plan->psi;
1735 nfft_init_guru(&(ths->set_nfft_plan_2d[r]),2,N,M,n,m, MALLOC_X| MALLOC_F| FFTW_INIT, FFTW_MEASURE);
1736 ths->set_nfft_plan_2d[r].flags = ths->set_nfft_plan_2d[r].flags | FG_PSI;
1737 ths->set_nfft_plan_2d[r].K=ths->act_nfft_plan->K;
1738 ths->set_nfft_plan_2d[r].psi=ths->act_nfft_plan->psi;
1739 }
1740 }
1741
1742 /* center plan */
1743 /* J/2 == floor(((double)J) / 2.0) */
1744 N[0]=X(exp2i)(J/2+1); n[0]=ths->sigma*N[0];
1745 N[1]=X(exp2i)(J/2+1); n[1]=ths->sigma*N[1];
1746 N[2]=X(exp2i)(J/2+1); n[2]=ths->sigma*N[2];
1747 nfft_init_guru(ths->center_nfft_plan,3,N,M,n, m, MALLOC_F| FFTW_INIT,
1748 FFTW_MEASURE);
1749 ths->center_nfft_plan->x= ths->act_nfft_plan->x;
1750 ths->center_nfft_plan->flags = ths->center_nfft_plan->flags|
1751 FG_PSI;
1752 ths->center_nfft_plan->K=ths->act_nfft_plan->K;
1753 ths->center_nfft_plan->psi=ths->act_nfft_plan->psi;
1754
1755 if(ths->flags & NSDFT)
1756 {
1757 ths->index_sparse_to_full=(int*)nfft_malloc(ths->N_total*sizeof(int));
1758 init_index_sparse_to_full_3d(ths);
1759 }
1760}
1761#endif
1762
1763#ifdef GAUSSIAN
1764void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
1765{
1766 ths->d=d;
1767
1768 if(ths->d==2)
1769 nsfft_init_2d(ths, J, M, m, flags);
1770 else
1771 nsfft_init_3d(ths, J, M, m, flags);
1772
1773
1774 ths->mv_trafo = (void (*) (void* ))nsfft_trafo;
1775 ths->mv_adjoint = (void (*) (void* ))nsfft_adjoint;
1776}
1777#else
1778void nsfft_init(nsfft_plan *ths, int d, int J, int M, int m, unsigned flags)
1779{
1780 UNUSED(ths);
1781 UNUSED(d);
1782 UNUSED(J);
1783 UNUSED(M);
1784 UNUSED(m);
1785 UNUSED(flags);
1786 fprintf(stderr,
1787 "\nError in kernel/nsfft_init: require GAUSSIAN window function\n");
1788}
1789#endif
1790
1791static void nsfft_finalize_2d(nsfft_plan *ths)
1792{
1793 int r;
1794
1795 if(ths->flags & NSDFT)
1797
1798 /* center plan */
1799 ths->center_nfft_plan->flags = ths->center_nfft_plan->flags ^ FG_PSI;
1800 nfft_finalize(ths->center_nfft_plan);
1801
1802 /* the 1d nffts */
1803 for(r=0;r<=X(log2i)(ths->act_nfft_plan->m);r++)
1804 {
1805 ths->set_nfft_plan_1d[r].flags =
1806 ths->set_nfft_plan_1d[r].flags ^ FG_PSI;
1807 nfft_finalize(&(ths->set_nfft_plan_1d[r]));
1808 }
1809
1810 /* finalize the small nffts */
1811 ths->act_nfft_plan->my_fftw_plan2=ths->set_fftw_plan2[0];
1812 ths->act_nfft_plan->my_fftw_plan1=ths->set_fftw_plan1[0];
1813
1814 for(r=1;r<=ths->J/2;r++)
1815 {
1816 fftw_destroy_plan(ths->set_fftw_plan2[r]);
1817 fftw_destroy_plan(ths->set_fftw_plan1[r]);
1818 }
1819
1820 /* r=0 */
1821 nfft_finalize(ths->act_nfft_plan);
1822
1824
1825 nfft_free(ths->set_fftw_plan2);
1826 nfft_free(ths->set_fftw_plan1);
1827
1828 nfft_free(ths->x_transposed);
1829
1830 nfft_free(ths->f_hat);
1831 nfft_free(ths->f);
1832}
1833
1834static void nsfft_finalize_3d(nsfft_plan *ths)
1835{
1836 int r;
1837
1838 if(ths->flags & NSDFT)
1840
1841 /* center plan */
1842 ths->center_nfft_plan->flags = ths->center_nfft_plan->flags ^ FG_PSI;
1843 nfft_finalize(ths->center_nfft_plan);
1844
1845 /* the 1d and 2d nffts */
1846 for(r=0;r<=X(log2i)(ths->act_nfft_plan->m);r++)
1847 {
1848 if(X(exp2i)(ths->J-r)>ths->act_nfft_plan->m)
1849 {
1850 ths->set_nfft_plan_2d[r].flags = ths->set_nfft_plan_2d[r].flags ^ FG_PSI;
1851 nfft_finalize(&(ths->set_nfft_plan_2d[r]));
1852 ths->set_nfft_plan_1d[r].flags = ths->set_nfft_plan_1d[r].flags ^ FG_PSI;
1853 nfft_finalize(&(ths->set_nfft_plan_1d[r]));
1854 }
1855 }
1856
1857 /* finalize the small nffts */
1858 ths->act_nfft_plan->my_fftw_plan2=ths->set_fftw_plan2[0];
1859 ths->act_nfft_plan->my_fftw_plan1=ths->set_fftw_plan1[0];
1860
1861 for(r=1;r<=(ths->J+1)/2;r++)
1862 {
1863 fftw_destroy_plan(ths->set_fftw_plan2[r]);
1864 fftw_destroy_plan(ths->set_fftw_plan1[r]);
1865 }
1866
1867 /* r=0 */
1868 nfft_finalize(ths->act_nfft_plan);
1869
1872
1873 nfft_free(ths->set_fftw_plan2);
1874 nfft_free(ths->set_fftw_plan1);
1875
1876 nfft_free(ths->x_102);
1877 nfft_free(ths->x_201);
1878 nfft_free(ths->x_120);
1879 nfft_free(ths->x_021);
1880
1881 nfft_free(ths->f_hat);
1882 nfft_free(ths->f);
1883}
1884
1885void nsfft_finalize(nsfft_plan *ths)
1886{
1887 if(ths->d==2)
1888 nsfft_finalize_2d(ths);
1889 else
1890 nsfft_finalize_3d(ths);
1891}
#define FG_PSI
Definition nfft3.h:188
#define MALLOC_X
Definition nfft3.h:193
#define PRE_ONE_PSI
Definition nfft3.h:200
#define MALLOC_F
Definition nfft3.h:195
#define FFTW_INIT
Definition nfft3.h:197
#define RSWAP(x, y)
Swap two vectors.
Definition infft.h:143
#define UNUSED(x)
Dummy use of unused parameters to silence compiler warnings.
Definition infft.h:1319
#define NSDFT
Definition nfft3.h:482
Internal header file for auxiliary definitions and functions.
Header file for the nfft3 library.
void nfft_vrand_unit_complex(fftw_complex *x, const NFFT_INT n)
Inits a vector of random complex numbers in .
void nfft_vrand_shifted_unit_double(double *x, const NFFT_INT n)
Inits a vector of random double numbers in .
void * nfft_malloc(size_t n)
void nfft_free(void *p)
data structure for an NSFFT (nonequispaced sparse fast Fourier transform) plan with double precision
Definition nfft3.h:478
nfft_plan * center_nfft_plan
central nfft block
Definition nfft3.h:478
nfft_plan * set_nfft_plan_2d
nfft plans for short nffts
Definition nfft3.h:478
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:478
int sigma
oversampling-factor
Definition nfft3.h:478
unsigned flags
flags for precomputation, malloc
Definition nfft3.h:478
fftw_complex * f
Samples.
Definition nfft3.h:478
int * index_sparse_to_full
index conversation, overflow for d=3, J=9!
Definition nfft3.h:478
nfft_plan * act_nfft_plan
current nfft block
Definition nfft3.h:478
void(* mv_adjoint)(void *)
Adjoint transform.
Definition nfft3.h:478
int d
dimension, rank; d = 2, 3
Definition nfft3.h:478
double * x_transposed
coordinate exchanged nodes, d = 2
Definition nfft3.h:478
nfft_plan * set_nfft_plan_1d
nfft plans for short nffts
Definition nfft3.h:478
int J
problem size, i.e., d=2: N_total=(J+4) 2^(J+1) d=3: N_total=2^J 6(2^((J+1)/2+1)-1)+2^(3(J/2+1))
Definition nfft3.h:478
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:478
NFFT_INT N_total
Total number of Fourier coefficients.
Definition nfft3.h:478
double * x_021
coordinate exchanged nodes, d=3
Definition nfft3.h:478
void(* mv_trafo)(void *)
Transform.
Definition nfft3.h:478