ITK  4.3.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 #include "itkMutexLockHolder.h"
27 
28 #if !defined(FFTW_WISDOM_ONLY)
29 // FFTW_WISDOM_ONLY is a "beyond guru" option that is only available in fftw 3.2.2
30 // to be compatible with all the fftw 3.x API, we need to define this away here:
31 #error "FFTW 3.3.2 or later is required so that FFTW_WISDOM_ONLY is defined."
32 #endif
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( 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  static PlanType Plan_dft_c2r_1d(int n,
72  ComplexType *in,
73  PixelType *out,
74  unsigned flags,
75  int threads=1,
76  bool canDestroyInput=false)
77  {
78  return Plan_dft_c2r(1, &n, in, out, flags, threads, canDestroyInput);
79  }
80 
81  static PlanType Plan_dft_c2r_2d(int nx,
82  int ny,
83  ComplexType *in,
84  PixelType *out,
85  unsigned flags,
86  int threads=1,
87  bool canDestroyInput=false)
88  {
89  int * sizes = new int[2];
90  sizes[0] = nx;
91  sizes[1] = ny;
92  PlanType plan = Plan_dft_c2r(2, sizes, in, out, flags, threads, canDestroyInput);
93  delete sizes;
94  return plan;
95  }
96 
97  static PlanType Plan_dft_c2r_3d(int nx,
98  int ny,
99  int nz,
100  ComplexType *in,
101  PixelType *out,
102  unsigned flags,
103  int threads=1,
104  bool canDestroyInput=false)
105  {
106  int * sizes = new int[3];
107  sizes[0] = nx;
108  sizes[1] = ny;
109  sizes[2] = nz;
110  PlanType plan = Plan_dft_c2r(3, sizes, in, out, flags, threads, canDestroyInput);
111  delete sizes;
112  return plan;
113  }
114 
115  static PlanType Plan_dft_c2r(int rank,
116  const int *n,
117  ComplexType *in,
118  PixelType *out,
119  unsigned flags,
120  int threads=1,
121  bool canDestroyInput=false)
122  {
123  MutexLockHolder< FFTWGlobalConfiguration::MutexType > lock( FFTWGlobalConfiguration::GetLockMutex() );
124  fftwf_plan_with_nthreads(threads);
125  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
126  // because FFTW_ESTIMATE guarantee to not destroy the input
127  unsigned roflags = flags;
128  if( ! (flags & FFTW_ESTIMATE) )
129  {
130  roflags = flags | FFTW_WISDOM_ONLY;
131  }
132  PlanType plan = fftwf_plan_dft_c2r(rank,n,in,out,roflags);
133  if( plan == NULL )
134  {
135  // no wisdom available for that plan
136  if( canDestroyInput )
137  {
138  // just create the plan
139  plan = fftwf_plan_dft_c2r(rank,n,in,out,flags);
140  }
141  else
142  {
143  // lets create a plan with a fake input to generate the wisdom
144  int total = 1;
145  for( int i=0; i<rank; i++ )
146  {
147  total *= n[i];
148  }
149  ComplexType * din = new ComplexType[total];
150  fftwf_plan_dft_c2r(rank,n,din,out,flags);
151  delete [] din;
152  // and then create the final plan - this time it shouldn't fail
153  plan = fftwf_plan_dft_c2r(rank,n,in,out,roflags);
154  }
155  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
156  }
157  assert( plan != NULL );
158  return plan;
159  }
160 
161 
162  static PlanType Plan_dft_r2c_1d(int n,
163  PixelType *in,
164  ComplexType *out,
165  unsigned flags,
166  int threads=1,
167  bool canDestroyInput=false)
168  {
169  return Plan_dft_r2c(1, &n, in, out, flags, threads, canDestroyInput);
170  }
171 
172  static PlanType Plan_dft_r2c_2d(int nx,
173  int ny,
174  PixelType *in,
175  ComplexType *out,
176  unsigned flags,
177  int threads=1,
178  bool canDestroyInput=false)
179  {
180  int * sizes = new int[2];
181  sizes[0] = nx;
182  sizes[1] = ny;
183  PlanType plan = Plan_dft_r2c(2, sizes, in, out, flags, threads, canDestroyInput);
184  delete sizes;
185  return plan;
186  }
187 
188  static PlanType Plan_dft_r2c_3d(int nx,
189  int ny,
190  int nz,
191  PixelType *in,
192  ComplexType *out,
193  unsigned flags,
194  int threads=1,
195  bool canDestroyInput=false)
196  {
197  int * sizes = new int[3];
198  sizes[0] = nx;
199  sizes[1] = ny;
200  sizes[2] = nz;
201  PlanType plan = Plan_dft_r2c(3, sizes, in, out, flags, threads, canDestroyInput);
202  delete sizes;
203  return plan;
204  }
205 
206  static PlanType Plan_dft_r2c(int rank,
207  const int *n,
208  PixelType *in,
209  ComplexType *out,
210  unsigned flags,
211  int threads=1,
212  bool canDestroyInput=false)
213  {
214  //
215  MutexLockHolder< FFTWGlobalConfiguration::MutexType > lock( FFTWGlobalConfiguration::GetLockMutex() );
216  fftwf_plan_with_nthreads(threads);
217  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
218  // because FFTW_ESTIMATE guarantee to not destroy the input
219  unsigned roflags = flags;
220  if( ! (flags & FFTW_ESTIMATE) )
221  {
222  roflags = flags | FFTW_WISDOM_ONLY;
223  }
224  PlanType plan = fftwf_plan_dft_r2c(rank,n,in,out,roflags);
225  if( plan == NULL )
226  {
227  // no wisdom available for that plan
228  if( canDestroyInput )
229  {
230  // just create the plan
231  plan = fftwf_plan_dft_r2c(rank,n,in,out,flags);
232  }
233  else
234  {
235  // lets create a plan with a fake input to generate the wisdom
236  int total = 1;
237  for( int i=0; i<rank; i++ )
238  {
239  total *= n[i];
240  }
241  PixelType * din = new PixelType[total];
242  fftwf_plan_dft_r2c(rank,n,din,out,flags);
243  delete [] din;
244  // and then create the final plan - this time it shouldn't fail
245  plan = fftwf_plan_dft_r2c(rank,n,in,out,roflags);
246  }
247  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
248  }
249  assert( plan != NULL );
250  return plan;
251  }
252 
253  static PlanType Plan_dft_1d(int n,
254  ComplexType *in,
255  ComplexType *out,
256  int sign,
257  unsigned flags,
258  int threads=1,
259  bool canDestroyInput=false)
260  {
261  return Plan_dft(1, &n, in, out,sign , flags, threads, canDestroyInput);
262  }
263 
264  static PlanType Plan_dft_2d(int nx,
265  int ny,
266  ComplexType *in,
267  ComplexType *out,
268  int sign,
269  unsigned flags,
270  int threads=1,
271  bool canDestroyInput=false)
272  {
273  int * sizes = new int[2];
274  sizes[0] = nx;
275  sizes[1] = ny;
276  PlanType plan = Plan_dft(2, sizes, in, out, sign, flags, threads, canDestroyInput);
277  delete sizes;
278  return plan;
279  }
280 
281  static PlanType Plan_dft_3d(int nx,
282  int ny,
283  int nz,
284  ComplexType *in,
285  ComplexType *out,
286  int sign,
287  unsigned flags,
288  int threads=1,
289  bool canDestroyInput=false)
290  {
291  int * sizes = new int[3];
292  sizes[0] = nx;
293  sizes[1] = ny;
294  sizes[2] = nz;
295  PlanType plan = Plan_dft(3, sizes, in, out, sign, flags, threads, canDestroyInput);
296  delete sizes;
297  return plan;
298  }
299 
300  static PlanType Plan_dft(int rank,
301  const int *n,
302  ComplexType *in,
303  ComplexType *out,
304  int sign,
305  unsigned flags,
306  int threads=1,
307  bool canDestroyInput=false)
308  {
309  MutexLockHolder< FFTWGlobalConfiguration::MutexType > lock( FFTWGlobalConfiguration::GetLockMutex() );
310  fftwf_plan_with_nthreads(threads);
311  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
312  // because FFTW_ESTIMATE guarantee to not destroy the input
313  unsigned roflags = flags;
314  if( ! (flags & FFTW_ESTIMATE) )
315  {
316  roflags = flags | FFTW_WISDOM_ONLY;
317  }
318  PlanType plan = fftwf_plan_dft(rank,n,in,out,sign,roflags);
319  if( plan == NULL )
320  {
321  // no wisdom available for that plan
322  if( canDestroyInput )
323  {
324  // just create the plan
325  plan = fftwf_plan_dft(rank,n,in,out,sign,flags);
326  }
327  else
328  {
329  // lets create a plan with a fake input to generate the wisdom
330  int total = 1;
331  for( int i=0; i<rank; i++ )
332  {
333  total *= n[i];
334  }
335  ComplexType * din = new ComplexType[total];
336  fftwf_plan_dft(rank,n,din,out,sign,flags);
337  delete [] din;
338  // and then create the final plan - this time it shouldn't fail
339  plan = fftwf_plan_dft(rank,n,in,out,sign,roflags);
340  }
341  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
342  }
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  MutexLockHolder< FFTWGlobalConfiguration::MutexType > lock( FFTWGlobalConfiguration::GetLockMutex() );
355  fftwf_destroy_plan(p);
356  }
357 };
358 
359 #endif // USE_FFTWF
360 
361 
362 #if defined( USE_FFTWD )
363 template< >
364 class Proxy< double >
365 {
366 public:
367  typedef double PixelType;
368  typedef fftw_complex ComplexType;
369  typedef fftw_plan PlanType;
370  typedef Proxy<double> Self;
371 
372  static PlanType Plan_dft_c2r_1d(int n,
373  ComplexType *in,
374  PixelType *out,
375  unsigned flags,
376  int threads=1,
377  bool canDestroyInput=false)
378  {
379  return Plan_dft_c2r(1, &n, in, out, flags, threads, canDestroyInput);
380  }
381 
382  static PlanType Plan_dft_c2r_2d(int nx,
383  int ny,
384  ComplexType *in,
385  PixelType *out,
386  unsigned flags,
387  int threads=1,
388  bool canDestroyInput=false)
389  {
390  int * sizes = new int[2];
391  sizes[0] = nx;
392  sizes[1] = ny;
393  PlanType plan = Plan_dft_c2r(2, sizes, in, out, flags, threads, canDestroyInput);
394  delete sizes;
395  return plan;
396  }
397 
398  static PlanType Plan_dft_c2r_3d(int nx,
399  int ny,
400  int nz,
401  ComplexType *in,
402  PixelType *out,
403  unsigned flags,
404  int threads=1,
405  bool canDestroyInput=false)
406  {
407  int * sizes = new int[3];
408  sizes[0] = nx;
409  sizes[1] = ny;
410  sizes[2] = nz;
411  PlanType plan = Plan_dft_c2r(3, sizes, in, out, flags, threads, canDestroyInput);
412  delete sizes;
413  return plan;
414  }
415 
416  static PlanType Plan_dft_c2r(int rank,
417  const int *n,
418  ComplexType *in,
419  PixelType *out,
420  unsigned flags,
421  int threads=1,
422  bool canDestroyInput=false)
423  {
424  MutexLockHolder< FFTWGlobalConfiguration::MutexType > lock( FFTWGlobalConfiguration::GetLockMutex() );
425  fftw_plan_with_nthreads(threads);
426  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
427  // because FFTW_ESTIMATE guarantee to not destroy the input
428  unsigned roflags = flags;
429  if( ! (flags & FFTW_ESTIMATE) )
430  {
431  roflags = flags | FFTW_WISDOM_ONLY;
432  }
433  PlanType plan = fftw_plan_dft_c2r(rank,n,in,out,roflags);
434  if( plan == NULL )
435  {
436  // no wisdom available for that plan
437  if( canDestroyInput )
438  {
439  // just create the plan
440  plan = fftw_plan_dft_c2r(rank,n,in,out,flags);
441  }
442  else
443  {
444  // lets create a plan with a fake input to generate the wisdom
445  int total = 1;
446  for( int i=0; i<rank; i++ )
447  {
448  total *= n[i];
449  }
450  ComplexType * din = new ComplexType[total];
451  fftw_plan_dft_c2r(rank,n,din,out,flags);
452  delete [] din;
453  // and then create the final plan - this time it shouldn't fail
454  plan = fftw_plan_dft_c2r(rank,n,in,out,roflags);
455  }
456  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
457  }
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  MutexLockHolder< FFTWGlobalConfiguration::MutexType > lock( FFTWGlobalConfiguration::GetLockMutex() );
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  assert( plan != NULL );
550  return plan;
551  }
552 
553  static PlanType Plan_dft_1d(int n,
554  ComplexType *in,
555  ComplexType *out,
556  int sign,
557  unsigned flags,
558  int threads=1,
559  bool canDestroyInput=false)
560  {
561  return Plan_dft(1, &n, in, out,sign , flags, threads, canDestroyInput);
562  }
563 
564  static PlanType Plan_dft_2d(int nx,
565  int ny,
566  ComplexType *in,
567  ComplexType *out,
568  int sign,
569  unsigned flags,
570  int threads=1,
571  bool canDestroyInput=false)
572  {
573  int * sizes = new int[2];
574  sizes[0] = nx;
575  sizes[1] = ny;
576  PlanType plan = Plan_dft(2, sizes, in, out, sign, flags, threads, canDestroyInput);
577  delete sizes;
578  return plan;
579  }
580 
581  static PlanType Plan_dft_3d(int nx,
582  int ny,
583  int nz,
584  ComplexType *in,
585  ComplexType *out,
586  int sign,
587  unsigned flags,
588  int threads=1,
589  bool canDestroyInput=false)
590  {
591  int * sizes = new int[3];
592  sizes[0] = nx;
593  sizes[1] = ny;
594  sizes[2] = nz;
595  PlanType plan = Plan_dft(3, sizes, in, out, sign, flags, threads, canDestroyInput);
596  delete sizes;
597  return plan;
598  }
599 
600  static PlanType Plan_dft(int rank,
601  const int *n,
602  ComplexType *in,
603  ComplexType *out,
604  int sign,
605  unsigned flags,
606  int threads=1,
607  bool canDestroyInput=false)
608  {
609  MutexLockHolder< FFTWGlobalConfiguration::MutexType > lock( FFTWGlobalConfiguration::GetLockMutex() );
610  fftw_plan_with_nthreads(threads);
611  // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE
612  // because FFTW_ESTIMATE guarantee to not destroy the input
613  unsigned roflags = flags;
614  if( ! (flags & FFTW_ESTIMATE) )
615  {
616  roflags = flags | FFTW_WISDOM_ONLY;
617  }
618  PlanType plan = fftw_plan_dft(rank,n,in,out,sign,roflags);
619  if( plan == NULL )
620  {
621  // no wisdom available for that plan
622  if( canDestroyInput )
623  {
624  // just create the plan
625  plan = fftw_plan_dft(rank,n,in,out,sign,flags);
626  }
627  else
628  {
629  // lets create a plan with a fake input to generate the wisdom
630  int total = 1;
631  for( int i=0; i<rank; i++ )
632  {
633  total *= n[i];
634  }
635  ComplexType * din = new ComplexType[total];
636  fftw_plan_dft(rank,n,din,out,sign,flags);
637  delete [] din;
638  // and then create the final plan - this time it shouldn't fail
639  plan = fftw_plan_dft(rank,n,in,out,sign,roflags);
640  }
641  FFTWGlobalConfiguration::SetNewWisdomAvailable(true);
642  }
643  assert( plan != NULL );
644  return plan;
645  }
646 
647 
648  static void Execute(PlanType p)
649  {
650  fftw_execute(p);
651  }
652  static void DestroyPlan(PlanType p)
653  {
654  MutexLockHolder< FFTWGlobalConfiguration::MutexType > lock( FFTWGlobalConfiguration::GetLockMutex() );
655  fftw_destroy_plan(p);
656  }
657 };
658 
659 #endif
660 }
661 }
662 #endif
663 
Tarballs of the nightly generated Doxygen documentation are available for the html, xml, and tag file.
Generated on Sun Dec 9 2012 00:59:43 for ITK by   doxygen 1.8.2