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