ITK  4.2.0
Insight Segmentation and Registration Toolkit
itkFFTWCommon.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef __itkFFTWCommon_h
19 #define __itkFFTWCommon_h
20 
21 #include "itkFFTWGlobalConfiguration.h"
22 #if defined( USE_FFTWF ) || defined( USE_FFTWD )
23 #include "fftw3.h"
24 #endif
25 
26 #if !defined(FFTW_WISDOM_ONLY)
27 // FFTW_WISDOM_ONLY is a "beyond guru" option that is only available in fftw 3.2.2
28 // to be compatible with all the fftw 3.x API, we need to define this away here:
29 #error "FFTW 3.3.2 or later is required so that FFTW_WISDOM_ONLY is defined."
30 #endif
31 
32 namespace itk
33 {
34 namespace fftw
35 {
48 template< typename TPixel >
49 class Proxy
50 {
51  // empty -- only double and float specializations work
52 protected:
53  Proxy() {};
54  ~Proxy() {};
55 };
57 
58 #if defined( USE_FFTWF )
59 
60 template< >
61 class Proxy< float >
62 {
63 public:
64  typedef float PixelType;
65  typedef fftwf_complex ComplexType;
66  typedef fftwf_plan PlanType;
67  typedef Proxy<float> Self;
68 
69  static PlanType Plan_dft_c2r_1d(int n,
70  ComplexType *in,
71  PixelType *out,
72  unsigned flags,
73  int threads=1,
74  bool canDestroyInput=false)
75  {
76  return Plan_dft_c2r(1, &n, in, out, flags, threads, canDestroyInput);
77  }
78 
79  static PlanType Plan_dft_c2r_2d(int nx,
80  int ny,
81  ComplexType *in,
82  PixelType *out,
83  unsigned flags,
84  int threads=1,
85  bool canDestroyInput=false)
86  {
87  int * sizes = new int[2];
88  sizes[0] = nx;
89  sizes[1] = ny;
90  PlanType plan = Plan_dft_c2r(2, sizes, in, out, flags, threads, canDestroyInput);
91  delete sizes;
92  return plan;
93  }
94 
95  static PlanType Plan_dft_c2r_3d(int nx,
96  int ny,
97  int nz,
98  ComplexType *in,
99  PixelType *out,
100  unsigned flags,
101  int threads=1,
102  bool canDestroyInput=false)
103  {
104  int * sizes = new int[3];
105  sizes[0] = nx;
106  sizes[1] = ny;
107  sizes[2] = nz;
108  PlanType plan = Plan_dft_c2r(3, sizes, in, out, flags, threads, canDestroyInput);
109  delete sizes;
110  return plan;
111  }
112 
113  static PlanType Plan_dft_c2r(int rank,
114  const int *n,
115  ComplexType *in,
116  PixelType *out,
117  unsigned flags,
118  int threads=1,
119  bool canDestroyInput=false)
120  {
121  FFTWGlobalConfiguration::Lock();
122  fftwf_plan_with_nthreads(threads);
123  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
124  // because FFTW_ESTIMATE guarantee to not destroy the input
125  unsigned roflags = flags;
126  if( ! (flags & FFTW_ESTIMATE) )
127  {
128  roflags = flags | FFTW_WISDOM_ONLY;
129  }
130  PlanType plan = fftwf_plan_dft_c2r(rank,n,in,out,roflags);
131  if( plan == NULL )
132  {
133  // no wisdom available for that plan
134  if( canDestroyInput )
135  {
136  // just create the plan
137  plan = fftwf_plan_dft_c2r(rank,n,in,out,flags);
138  }
139  else
140  {
141  // lets create a plan with a fake input to generate the wisdom
142  int total = 1;
143  for( int i=0; i<rank; i++ )
144  {
145  total *= n[i];
146  }
147  ComplexType * din = new ComplexType[total];
148  fftwf_plan_dft_c2r(rank,n,din,out,flags);
149  delete [] din;
150  // and then create the final plan - this time it shouldn't fail
151  plan = fftwf_plan_dft_c2r(rank,n,in,out,roflags);
152  }
153  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
154  }
155  FFTWGlobalConfiguration::Unlock();
156  assert( plan != NULL );
157  return plan;
158  }
159 
160 
161  static PlanType Plan_dft_r2c_1d(int n,
162  PixelType *in,
163  ComplexType *out,
164  unsigned flags,
165  int threads=1,
166  bool canDestroyInput=false)
167  {
168  return Plan_dft_r2c(1, &n, in, out, flags, threads, canDestroyInput);
169  }
170 
171  static PlanType Plan_dft_r2c_2d(int nx,
172  int ny,
173  PixelType *in,
174  ComplexType *out,
175  unsigned flags,
176  int threads=1,
177  bool canDestroyInput=false)
178  {
179  int * sizes = new int[2];
180  sizes[0] = nx;
181  sizes[1] = ny;
182  PlanType plan = Plan_dft_r2c(2, sizes, in, out, flags, threads, canDestroyInput);
183  delete sizes;
184  return plan;
185  }
186 
187  static PlanType Plan_dft_r2c_3d(int nx,
188  int ny,
189  int nz,
190  PixelType *in,
191  ComplexType *out,
192  unsigned flags,
193  int threads=1,
194  bool canDestroyInput=false)
195  {
196  int * sizes = new int[3];
197  sizes[0] = nx;
198  sizes[1] = ny;
199  sizes[2] = nz;
200  PlanType plan = Plan_dft_r2c(3, sizes, in, out, flags, threads, canDestroyInput);
201  delete sizes;
202  return plan;
203  }
204 
205  static PlanType Plan_dft_r2c(int rank,
206  const int *n,
207  PixelType *in,
208  ComplexType *out,
209  unsigned flags,
210  int threads=1,
211  bool canDestroyInput=false)
212  {
213  FFTWGlobalConfiguration::Lock();
214  fftwf_plan_with_nthreads(threads);
215  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
216  // because FFTW_ESTIMATE guarantee to not destroy the input
217  unsigned roflags = flags;
218  if( ! (flags & FFTW_ESTIMATE) )
219  {
220  roflags = flags | FFTW_WISDOM_ONLY;
221  }
222  PlanType plan = fftwf_plan_dft_r2c(rank,n,in,out,roflags);
223  if( plan == NULL )
224  {
225  // no wisdom available for that plan
226  if( canDestroyInput )
227  {
228  // just create the plan
229  plan = fftwf_plan_dft_r2c(rank,n,in,out,flags);
230  }
231  else
232  {
233  // lets create a plan with a fake input to generate the wisdom
234  int total = 1;
235  for( int i=0; i<rank; i++ )
236  {
237  total *= n[i];
238  }
239  PixelType * din = new PixelType[total];
240  fftwf_plan_dft_r2c(rank,n,din,out,flags);
241  delete [] din;
242  // and then create the final plan - this time it shouldn't fail
243  plan = fftwf_plan_dft_r2c(rank,n,in,out,roflags);
244  }
245  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
246  }
247  FFTWGlobalConfiguration::Unlock();
248  assert( plan != NULL );
249  return plan;
250  }
251 
252  static PlanType Plan_dft_1d(int n,
253  ComplexType *in,
254  ComplexType *out,
255  int sign,
256  unsigned flags,
257  int threads=1,
258  bool canDestroyInput=false)
259  {
260  return Plan_dft(1, &n, in, out,sign , flags, threads, canDestroyInput);
261  }
262 
263  static PlanType Plan_dft_2d(int nx,
264  int ny,
265  ComplexType *in,
266  ComplexType *out,
267  int sign,
268  unsigned flags,
269  int threads=1,
270  bool canDestroyInput=false)
271  {
272  int * sizes = new int[2];
273  sizes[0] = nx;
274  sizes[1] = ny;
275  PlanType plan = Plan_dft(2, sizes, in, out, sign, flags, threads, canDestroyInput);
276  delete sizes;
277  return plan;
278  }
279 
280  static PlanType Plan_dft_3d(int nx,
281  int ny,
282  int nz,
283  ComplexType *in,
284  ComplexType *out,
285  int sign,
286  unsigned flags,
287  int threads=1,
288  bool canDestroyInput=false)
289  {
290  int * sizes = new int[3];
291  sizes[0] = nx;
292  sizes[1] = ny;
293  sizes[2] = nz;
294  PlanType plan = Plan_dft(3, sizes, in, out, sign, flags, threads, canDestroyInput);
295  delete sizes;
296  return plan;
297  }
298 
299  static PlanType Plan_dft(int rank,
300  const int *n,
301  ComplexType *in,
302  ComplexType *out,
303  int sign,
304  unsigned flags,
305  int threads=1,
306  bool canDestroyInput=false)
307  {
308  FFTWGlobalConfiguration::Lock();
309  fftwf_plan_with_nthreads(threads);
310  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
311  // because FFTW_ESTIMATE guarantee to not destroy the input
312  unsigned roflags = flags;
313  if( ! (flags & FFTW_ESTIMATE) )
314  {
315  roflags = flags | FFTW_WISDOM_ONLY;
316  }
317  PlanType plan = fftwf_plan_dft(rank,n,in,out,sign,roflags);
318  if( plan == NULL )
319  {
320  // no wisdom available for that plan
321  if( canDestroyInput )
322  {
323  // just create the plan
324  plan = fftwf_plan_dft(rank,n,in,out,sign,flags);
325  }
326  else
327  {
328  // lets create a plan with a fake input to generate the wisdom
329  int total = 1;
330  for( int i=0; i<rank; i++ )
331  {
332  total *= n[i];
333  }
334  ComplexType * din = new ComplexType[total];
335  fftwf_plan_dft(rank,n,din,out,sign,flags);
336  delete [] din;
337  // and then create the final plan - this time it shouldn't fail
338  plan = fftwf_plan_dft(rank,n,in,out,sign,roflags);
339  }
340  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
341  }
342  FFTWGlobalConfiguration::Unlock();
343  assert( plan != NULL );
344  return plan;
345  }
346 
347 
348  static void Execute(PlanType p)
349  {
350  fftwf_execute(p);
351  }
352  static void DestroyPlan(PlanType p)
353  {
354  fftwf_destroy_plan(p);
355  }
356 };
357 
358 #endif // USE_FFTWF
359 
360 
361 #if defined( USE_FFTWD )
362 template< >
363 class Proxy< double >
364 {
365 public:
366  typedef double PixelType;
367  typedef fftw_complex ComplexType;
368  typedef fftw_plan PlanType;
369  typedef Proxy<double> Self;
370 
371  static PlanType Plan_dft_c2r_1d(int n,
372  ComplexType *in,
373  PixelType *out,
374  unsigned flags,
375  int threads=1,
376  bool canDestroyInput=false)
377  {
378  return Plan_dft_c2r(1, &n, in, out, flags, threads, canDestroyInput);
379  }
380 
381  static PlanType Plan_dft_c2r_2d(int nx,
382  int ny,
383  ComplexType *in,
384  PixelType *out,
385  unsigned flags,
386  int threads=1,
387  bool canDestroyInput=false)
388  {
389  int * sizes = new int[2];
390  sizes[0] = nx;
391  sizes[1] = ny;
392  PlanType plan = Plan_dft_c2r(2, sizes, in, out, flags, threads, canDestroyInput);
393  delete sizes;
394  return plan;
395  }
396 
397  static PlanType Plan_dft_c2r_3d(int nx,
398  int ny,
399  int nz,
400  ComplexType *in,
401  PixelType *out,
402  unsigned flags,
403  int threads=1,
404  bool canDestroyInput=false)
405  {
406  int * sizes = new int[3];
407  sizes[0] = nx;
408  sizes[1] = ny;
409  sizes[2] = nz;
410  PlanType plan = Plan_dft_c2r(3, sizes, in, out, flags, threads, canDestroyInput);
411  delete sizes;
412  return plan;
413  }
414 
415  static PlanType Plan_dft_c2r(int rank,
416  const int *n,
417  ComplexType *in,
418  PixelType *out,
419  unsigned flags,
420  int threads=1,
421  bool canDestroyInput=false)
422  {
423  FFTWGlobalConfiguration::Lock();
424  fftw_plan_with_nthreads(threads);
425  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
426  // because FFTW_ESTIMATE guarantee to not destroy the input
427  unsigned roflags = flags;
428  if( ! (flags & FFTW_ESTIMATE) )
429  {
430  roflags = flags | FFTW_WISDOM_ONLY;
431  }
432  PlanType plan = fftw_plan_dft_c2r(rank,n,in,out,roflags);
433  if( plan == NULL )
434  {
435  // no wisdom available for that plan
436  if( canDestroyInput )
437  {
438  // just create the plan
439  plan = fftw_plan_dft_c2r(rank,n,in,out,flags);
440  }
441  else
442  {
443  // lets create a plan with a fake input to generate the wisdom
444  int total = 1;
445  for( int i=0; i<rank; i++ )
446  {
447  total *= n[i];
448  }
449  ComplexType * din = new ComplexType[total];
450  fftw_plan_dft_c2r(rank,n,din,out,flags);
451  delete [] din;
452  // and then create the final plan - this time it shouldn't fail
453  plan = fftw_plan_dft_c2r(rank,n,in,out,roflags);
454  }
455  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
456  }
457  FFTWGlobalConfiguration::Unlock();
458  assert( plan != NULL );
459  return plan;
460  }
461 
462 
463  static PlanType Plan_dft_r2c_1d(int n,
464  PixelType *in,
465  ComplexType *out,
466  unsigned flags,
467  int threads=1,
468  bool canDestroyInput=false)
469  {
470  return Plan_dft_r2c(1, &n, in, out, flags, threads, canDestroyInput);
471  }
472 
473  static PlanType Plan_dft_r2c_2d(int nx,
474  int ny,
475  PixelType *in,
476  ComplexType *out,
477  unsigned flags,
478  int threads=1,
479  bool canDestroyInput=false)
480  {
481  int * sizes = new int[2];
482  sizes[0] = nx;
483  sizes[1] = ny;
484  PlanType plan = Plan_dft_r2c(2, sizes, in, out, flags, threads, canDestroyInput);
485  delete sizes;
486  return plan;
487  }
488 
489  static PlanType Plan_dft_r2c_3d(int nx,
490  int ny,
491  int nz,
492  PixelType *in,
493  ComplexType *out,
494  unsigned flags,
495  int threads=1,
496  bool canDestroyInput=false)
497  {
498  int * sizes = new int[3];
499  sizes[0] = nx;
500  sizes[1] = ny;
501  sizes[2] = nz;
502  PlanType plan = Plan_dft_r2c(3, sizes, in, out, flags, threads, canDestroyInput);
503  delete sizes;
504  return plan;
505  }
506 
507  static PlanType Plan_dft_r2c(int rank,
508  const int *n,
509  PixelType *in,
510  ComplexType *out,
511  unsigned flags,
512  int threads=1,
513  bool canDestroyInput=false)
514  {
515  FFTWGlobalConfiguration::Lock();
516  fftw_plan_with_nthreads(threads);
517  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
518  // because FFTW_ESTIMATE guarantee to not destroy the input
519  unsigned roflags = flags;
520  if( ! (flags & FFTW_ESTIMATE) )
521  {
522  roflags = flags | FFTW_WISDOM_ONLY;
523  }
524  PlanType plan = fftw_plan_dft_r2c(rank,n,in,out,roflags);
525  if( plan == NULL )
526  {
527  // no wisdom available for that plan
528  if( canDestroyInput )
529  {
530  // just create the plan
531  plan = fftw_plan_dft_r2c(rank,n,in,out,flags);
532  }
533  else
534  {
535  // lets create a plan with a fake input to generate the wisdom
536  int total = 1;
537  for( int i=0; i<rank; i++ )
538  {
539  total *= n[i];
540  }
541  PixelType * din = new PixelType[total];
542  fftw_plan_dft_r2c(rank,n,din,out,flags);
543  delete [] din;
544  // and then create the final plan - this time it shouldn't fail
545  plan = fftw_plan_dft_r2c(rank,n,in,out,roflags);
546  }
547  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
548  }
549  FFTWGlobalConfiguration::Unlock();
550  assert( plan != NULL );
551  return plan;
552  }
553 
554  static PlanType Plan_dft_1d(int n,
555  ComplexType *in,
556  ComplexType *out,
557  int sign,
558  unsigned flags,
559  int threads=1,
560  bool canDestroyInput=false)
561  {
562  return Plan_dft(1, &n, in, out,sign , flags, threads, canDestroyInput);
563  }
564 
565  static PlanType Plan_dft_2d(int nx,
566  int ny,
567  ComplexType *in,
568  ComplexType *out,
569  int sign,
570  unsigned flags,
571  int threads=1,
572  bool canDestroyInput=false)
573  {
574  int * sizes = new int[2];
575  sizes[0] = nx;
576  sizes[1] = ny;
577  PlanType plan = Plan_dft(2, sizes, in, out, sign, flags, threads, canDestroyInput);
578  delete sizes;
579  return plan;
580  }
581 
582  static PlanType Plan_dft_3d(int nx,
583  int ny,
584  int nz,
585  ComplexType *in,
586  ComplexType *out,
587  int sign,
588  unsigned flags,
589  int threads=1,
590  bool canDestroyInput=false)
591  {
592  int * sizes = new int[3];
593  sizes[0] = nx;
594  sizes[1] = ny;
595  sizes[2] = nz;
596  PlanType plan = Plan_dft(3, sizes, in, out, sign, flags, threads, canDestroyInput);
597  delete sizes;
598  return plan;
599  }
600 
601  static PlanType Plan_dft(int rank,
602  const int *n,
603  ComplexType *in,
604  ComplexType *out,
605  int sign,
606  unsigned flags,
607  int threads=1,
608  bool canDestroyInput=false)
609  {
610  FFTWGlobalConfiguration::Lock();
611  fftw_plan_with_nthreads(threads);
612  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
613  // because FFTW_ESTIMATE guarantee to not destroy the input
614  unsigned roflags = flags;
615  if( ! (flags & FFTW_ESTIMATE) )
616  {
617  roflags = flags | FFTW_WISDOM_ONLY;
618  }
619  PlanType plan = fftw_plan_dft(rank,n,in,out,sign,roflags);
620  if( plan == NULL )
621  {
622  // no wisdom available for that plan
623  if( canDestroyInput )
624  {
625  // just create the plan
626  plan = fftw_plan_dft(rank,n,in,out,sign,flags);
627  }
628  else
629  {
630  // lets create a plan with a fake input to generate the wisdom
631  int total = 1;
632  for( int i=0; i<rank; i++ )
633  {
634  total *= n[i];
635  }
636  ComplexType * din = new ComplexType[total];
637  fftw_plan_dft(rank,n,din,out,sign,flags);
638  delete [] din;
639  // and then create the final plan - this time it shouldn't fail
640  plan = fftw_plan_dft(rank,n,in,out,sign,roflags);
641  }
642  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
643  }
644  FFTWGlobalConfiguration::Unlock();
645  assert( plan != NULL );
646  return plan;
647  }
648 
649 
650  static void Execute(PlanType p)
651  {
652  fftw_execute(p);
653  }
654  static void DestroyPlan(PlanType p)
655  {
656  fftw_destroy_plan(p);
657  }
658 };
659 
660 #endif
661 }
662 }
663 #endif
664 
Tarballs of the nightly generated Doxygen documentation are available for the html, xml, and tag file.
Generated on Tue Jul 10 2012 23:25:25 for ITK by doxygen 1.8.1