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