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