 |
ITK
4.1.0
Insight Segmentation and Registration Toolkit
|
|
ITK
4.1.0
Insight Segmentation and Registration Toolkit
|
00001 /*========================================================================= 00002 * 00003 * Copyright Insight Software Consortium 00004 * 00005 * Licensed under the Apache License, Version 2.0 (the "License"); 00006 * you may not use this file except in compliance with the License. 00007 * You may obtain a copy of the License at 00008 * 00009 * http://www.apache.org/licenses/LICENSE-2.0.txt 00010 * 00011 * Unless required by applicable law or agreed to in writing, software 00012 * distributed under the License is distributed on an "AS IS" BASIS, 00013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 00014 * See the License for the specific language governing permissions and 00015 * limitations under the License. 00016 * 00017 *=========================================================================*/ 00018 #ifndef __itkFFTWCommon_h 00019 #define __itkFFTWCommon_h 00020 00021 #include "itkFFTWGlobalConfiguration.h" 00022 #if defined( USE_FFTWF ) || defined( USE_FFTWD ) 00023 #include "fftw3.h" 00024 #endif 00025 00026 #if !defined(FFTW_WISDOM_ONLY) 00027 // FFTW_WISDOM_ONLY is a "beyond guru" option that is only available in fftw 3.2.2 00028 // to be compatible with all the fftw 3.x API, we need to define this away here: 00029 #error "FFTW 3.3.2 or later is required so that FFTW_WISDOM_ONLY is defined." 00030 #endif 00031 00032 namespace itk 00033 { 00034 namespace fftw 00035 { 00048 template< typename TPixel > 00049 class Proxy 00050 { 00051 // empty -- only double and float specializations work 00052 protected: 00053 Proxy() {}; 00054 ~Proxy() {}; 00055 }; 00057 00058 #if defined( USE_FFTWF ) 00059 00060 template< > 00061 class Proxy< float > 00062 { 00063 public: 00064 typedef float PixelType; 00065 typedef fftwf_complex ComplexType; 00066 typedef fftwf_plan PlanType; 00067 typedef Proxy<float> Self; 00068 00069 static PlanType Plan_dft_c2r_1d(int n, 00070 ComplexType *in, 00071 PixelType *out, 00072 unsigned flags, 00073 int threads=1, 00074 bool canDestroyInput=false) 00075 { 00076 return Plan_dft_c2r(1, &n, in, out, flags, threads, canDestroyInput); 00077 } 00078 00079 static PlanType Plan_dft_c2r_2d(int nx, 00080 int ny, 00081 ComplexType *in, 00082 PixelType *out, 00083 unsigned flags, 00084 int threads=1, 00085 bool canDestroyInput=false) 00086 { 00087 int * sizes = new int[2]; 00088 sizes[0] = nx; 00089 sizes[1] = ny; 00090 PlanType plan = Plan_dft_c2r(2, sizes, in, out, flags, threads, canDestroyInput); 00091 delete sizes; 00092 return plan; 00093 } 00094 00095 static PlanType Plan_dft_c2r_3d(int nx, 00096 int ny, 00097 int nz, 00098 ComplexType *in, 00099 PixelType *out, 00100 unsigned flags, 00101 int threads=1, 00102 bool canDestroyInput=false) 00103 { 00104 int * sizes = new int[3]; 00105 sizes[0] = nx; 00106 sizes[1] = ny; 00107 sizes[2] = nz; 00108 PlanType plan = Plan_dft_c2r(3, sizes, in, out, flags, threads, canDestroyInput); 00109 delete sizes; 00110 return plan; 00111 } 00112 00113 static PlanType Plan_dft_c2r(int rank, 00114 const int *n, 00115 ComplexType *in, 00116 PixelType *out, 00117 unsigned flags, 00118 int threads=1, 00119 bool canDestroyInput=false) 00120 { 00121 FFTWGlobalConfiguration::Lock(); 00122 fftwf_plan_with_nthreads(threads); 00123 // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE 00124 // because FFTW_ESTIMATE guarantee to not destroy the input 00125 unsigned roflags = flags; 00126 if( ! (flags & FFTW_ESTIMATE) ) 00127 { 00128 roflags = flags | FFTW_WISDOM_ONLY; 00129 } 00130 PlanType plan = fftwf_plan_dft_c2r(rank,n,in,out,roflags); 00131 if( plan == NULL ) 00132 { 00133 // no wisdom available for that plan 00134 if( canDestroyInput ) 00135 { 00136 // just create the plan 00137 plan = fftwf_plan_dft_c2r(rank,n,in,out,flags); 00138 } 00139 else 00140 { 00141 // lets create a plan with a fake input to generate the wisdom 00142 int total = 1; 00143 for( int i=0; i<rank; i++ ) 00144 { 00145 total *= n[i]; 00146 } 00147 ComplexType * din = new ComplexType[total]; 00148 fftwf_plan_dft_c2r(rank,n,din,out,flags); 00149 delete [] din; 00150 // and then create the final plan - this time it shouldn't fail 00151 plan = fftwf_plan_dft_c2r(rank,n,in,out,roflags); 00152 } 00153 FFTWGlobalConfiguration::SetNewWisdomAvailable(true); 00154 } 00155 FFTWGlobalConfiguration::Unlock(); 00156 assert( plan != NULL ); 00157 return plan; 00158 } 00159 00160 00161 static PlanType Plan_dft_r2c_1d(int n, 00162 PixelType *in, 00163 ComplexType *out, 00164 unsigned flags, 00165 int threads=1, 00166 bool canDestroyInput=false) 00167 { 00168 return Plan_dft_r2c(1, &n, in, out, flags, threads, canDestroyInput); 00169 } 00170 00171 static PlanType Plan_dft_r2c_2d(int nx, 00172 int ny, 00173 PixelType *in, 00174 ComplexType *out, 00175 unsigned flags, 00176 int threads=1, 00177 bool canDestroyInput=false) 00178 { 00179 int * sizes = new int[2]; 00180 sizes[0] = nx; 00181 sizes[1] = ny; 00182 PlanType plan = Plan_dft_r2c(2, sizes, in, out, flags, threads, canDestroyInput); 00183 delete sizes; 00184 return plan; 00185 } 00186 00187 static PlanType Plan_dft_r2c_3d(int nx, 00188 int ny, 00189 int nz, 00190 PixelType *in, 00191 ComplexType *out, 00192 unsigned flags, 00193 int threads=1, 00194 bool canDestroyInput=false) 00195 { 00196 int * sizes = new int[3]; 00197 sizes[0] = nx; 00198 sizes[1] = ny; 00199 sizes[2] = nz; 00200 PlanType plan = Plan_dft_r2c(3, sizes, in, out, flags, threads, canDestroyInput); 00201 delete sizes; 00202 return plan; 00203 } 00204 00205 static PlanType Plan_dft_r2c(int rank, 00206 const int *n, 00207 PixelType *in, 00208 ComplexType *out, 00209 unsigned flags, 00210 int threads=1, 00211 bool canDestroyInput=false) 00212 { 00213 FFTWGlobalConfiguration::Lock(); 00214 fftwf_plan_with_nthreads(threads); 00215 // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE 00216 // because FFTW_ESTIMATE guarantee to not destroy the input 00217 unsigned roflags = flags; 00218 if( ! (flags & FFTW_ESTIMATE) ) 00219 { 00220 roflags = flags | FFTW_WISDOM_ONLY; 00221 } 00222 PlanType plan = fftwf_plan_dft_r2c(rank,n,in,out,roflags); 00223 if( plan == NULL ) 00224 { 00225 // no wisdom available for that plan 00226 if( canDestroyInput ) 00227 { 00228 // just create the plan 00229 plan = fftwf_plan_dft_r2c(rank,n,in,out,flags); 00230 } 00231 else 00232 { 00233 // lets create a plan with a fake input to generate the wisdom 00234 int total = 1; 00235 for( int i=0; i<rank; i++ ) 00236 { 00237 total *= n[i]; 00238 } 00239 PixelType * din = new PixelType[total]; 00240 fftwf_plan_dft_r2c(rank,n,din,out,flags); 00241 delete [] din; 00242 // and then create the final plan - this time it shouldn't fail 00243 plan = fftwf_plan_dft_r2c(rank,n,in,out,roflags); 00244 } 00245 FFTWGlobalConfiguration::SetNewWisdomAvailable(true); 00246 } 00247 FFTWGlobalConfiguration::Unlock(); 00248 assert( plan != NULL ); 00249 return plan; 00250 } 00251 00252 static PlanType Plan_dft_1d(int n, 00253 ComplexType *in, 00254 ComplexType *out, 00255 int sign, 00256 unsigned flags, 00257 int threads=1, 00258 bool canDestroyInput=false) 00259 { 00260 return Plan_dft(1, &n, in, out,sign , flags, threads, canDestroyInput); 00261 } 00262 00263 static PlanType Plan_dft_2d(int nx, 00264 int ny, 00265 ComplexType *in, 00266 ComplexType *out, 00267 int sign, 00268 unsigned flags, 00269 int threads=1, 00270 bool canDestroyInput=false) 00271 { 00272 int * sizes = new int[2]; 00273 sizes[0] = nx; 00274 sizes[1] = ny; 00275 PlanType plan = Plan_dft(2, sizes, in, out, sign, flags, threads, canDestroyInput); 00276 delete sizes; 00277 return plan; 00278 } 00279 00280 static PlanType Plan_dft_3d(int nx, 00281 int ny, 00282 int nz, 00283 ComplexType *in, 00284 ComplexType *out, 00285 int sign, 00286 unsigned flags, 00287 int threads=1, 00288 bool canDestroyInput=false) 00289 { 00290 int * sizes = new int[3]; 00291 sizes[0] = nx; 00292 sizes[1] = ny; 00293 sizes[2] = nz; 00294 PlanType plan = Plan_dft(3, sizes, in, out, sign, flags, threads, canDestroyInput); 00295 delete sizes; 00296 return plan; 00297 } 00298 00299 static PlanType Plan_dft(int rank, 00300 const int *n, 00301 ComplexType *in, 00302 ComplexType *out, 00303 int sign, 00304 unsigned flags, 00305 int threads=1, 00306 bool canDestroyInput=false) 00307 { 00308 FFTWGlobalConfiguration::Lock(); 00309 fftwf_plan_with_nthreads(threads); 00310 // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE 00311 // because FFTW_ESTIMATE guarantee to not destroy the input 00312 unsigned roflags = flags; 00313 if( ! (flags & FFTW_ESTIMATE) ) 00314 { 00315 roflags = flags | FFTW_WISDOM_ONLY; 00316 } 00317 PlanType plan = fftwf_plan_dft(rank,n,in,out,sign,roflags); 00318 if( plan == NULL ) 00319 { 00320 // no wisdom available for that plan 00321 if( canDestroyInput ) 00322 { 00323 // just create the plan 00324 plan = fftwf_plan_dft(rank,n,in,out,sign,flags); 00325 } 00326 else 00327 { 00328 // lets create a plan with a fake input to generate the wisdom 00329 int total = 1; 00330 for( int i=0; i<rank; i++ ) 00331 { 00332 total *= n[i]; 00333 } 00334 ComplexType * din = new ComplexType[total]; 00335 fftwf_plan_dft(rank,n,din,out,sign,flags); 00336 delete [] din; 00337 // and then create the final plan - this time it shouldn't fail 00338 plan = fftwf_plan_dft(rank,n,in,out,sign,roflags); 00339 } 00340 FFTWGlobalConfiguration::SetNewWisdomAvailable(true); 00341 } 00342 FFTWGlobalConfiguration::Unlock(); 00343 assert( plan != NULL ); 00344 return plan; 00345 } 00346 00347 00348 static void Execute(PlanType p) 00349 { 00350 fftwf_execute(p); 00351 } 00352 static void DestroyPlan(PlanType p) 00353 { 00354 fftwf_destroy_plan(p); 00355 } 00356 }; 00357 00358 #endif // USE_FFTWF 00359 00360 00361 #if defined( USE_FFTWD ) 00362 template< > 00363 class Proxy< double > 00364 { 00365 public: 00366 typedef double PixelType; 00367 typedef fftw_complex ComplexType; 00368 typedef fftw_plan PlanType; 00369 typedef Proxy<double> Self; 00370 00371 static PlanType Plan_dft_c2r_1d(int n, 00372 ComplexType *in, 00373 PixelType *out, 00374 unsigned flags, 00375 int threads=1, 00376 bool canDestroyInput=false) 00377 { 00378 return Plan_dft_c2r(1, &n, in, out, flags, threads, canDestroyInput); 00379 } 00380 00381 static PlanType Plan_dft_c2r_2d(int nx, 00382 int ny, 00383 ComplexType *in, 00384 PixelType *out, 00385 unsigned flags, 00386 int threads=1, 00387 bool canDestroyInput=false) 00388 { 00389 int * sizes = new int[2]; 00390 sizes[0] = nx; 00391 sizes[1] = ny; 00392 PlanType plan = Plan_dft_c2r(2, sizes, in, out, flags, threads, canDestroyInput); 00393 delete sizes; 00394 return plan; 00395 } 00396 00397 static PlanType Plan_dft_c2r_3d(int nx, 00398 int ny, 00399 int nz, 00400 ComplexType *in, 00401 PixelType *out, 00402 unsigned flags, 00403 int threads=1, 00404 bool canDestroyInput=false) 00405 { 00406 int * sizes = new int[3]; 00407 sizes[0] = nx; 00408 sizes[1] = ny; 00409 sizes[2] = nz; 00410 PlanType plan = Plan_dft_c2r(3, sizes, in, out, flags, threads, canDestroyInput); 00411 delete sizes; 00412 return plan; 00413 } 00414 00415 static PlanType Plan_dft_c2r(int rank, 00416 const int *n, 00417 ComplexType *in, 00418 PixelType *out, 00419 unsigned flags, 00420 int threads=1, 00421 bool canDestroyInput=false) 00422 { 00423 FFTWGlobalConfiguration::Lock(); 00424 fftw_plan_with_nthreads(threads); 00425 // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE 00426 // because FFTW_ESTIMATE guarantee to not destroy the input 00427 unsigned roflags = flags; 00428 if( ! (flags & FFTW_ESTIMATE) ) 00429 { 00430 roflags = flags | FFTW_WISDOM_ONLY; 00431 } 00432 PlanType plan = fftw_plan_dft_c2r(rank,n,in,out,roflags); 00433 if( plan == NULL ) 00434 { 00435 // no wisdom available for that plan 00436 if( canDestroyInput ) 00437 { 00438 // just create the plan 00439 plan = fftw_plan_dft_c2r(rank,n,in,out,flags); 00440 } 00441 else 00442 { 00443 // lets create a plan with a fake input to generate the wisdom 00444 int total = 1; 00445 for( int i=0; i<rank; i++ ) 00446 { 00447 total *= n[i]; 00448 } 00449 ComplexType * din = new ComplexType[total]; 00450 fftw_plan_dft_c2r(rank,n,din,out,flags); 00451 delete [] din; 00452 // and then create the final plan - this time it shouldn't fail 00453 plan = fftw_plan_dft_c2r(rank,n,in,out,roflags); 00454 } 00455 FFTWGlobalConfiguration::SetNewWisdomAvailable(true); 00456 } 00457 FFTWGlobalConfiguration::Unlock(); 00458 assert( plan != NULL ); 00459 return plan; 00460 } 00461 00462 00463 static PlanType Plan_dft_r2c_1d(int n, 00464 PixelType *in, 00465 ComplexType *out, 00466 unsigned flags, 00467 int threads=1, 00468 bool canDestroyInput=false) 00469 { 00470 return Plan_dft_r2c(1, &n, in, out, flags, threads, canDestroyInput); 00471 } 00472 00473 static PlanType Plan_dft_r2c_2d(int nx, 00474 int ny, 00475 PixelType *in, 00476 ComplexType *out, 00477 unsigned flags, 00478 int threads=1, 00479 bool canDestroyInput=false) 00480 { 00481 int * sizes = new int[2]; 00482 sizes[0] = nx; 00483 sizes[1] = ny; 00484 PlanType plan = Plan_dft_r2c(2, sizes, in, out, flags, threads, canDestroyInput); 00485 delete sizes; 00486 return plan; 00487 } 00488 00489 static PlanType Plan_dft_r2c_3d(int nx, 00490 int ny, 00491 int nz, 00492 PixelType *in, 00493 ComplexType *out, 00494 unsigned flags, 00495 int threads=1, 00496 bool canDestroyInput=false) 00497 { 00498 int * sizes = new int[3]; 00499 sizes[0] = nx; 00500 sizes[1] = ny; 00501 sizes[2] = nz; 00502 PlanType plan = Plan_dft_r2c(3, sizes, in, out, flags, threads, canDestroyInput); 00503 delete sizes; 00504 return plan; 00505 } 00506 00507 static PlanType Plan_dft_r2c(int rank, 00508 const int *n, 00509 PixelType *in, 00510 ComplexType *out, 00511 unsigned flags, 00512 int threads=1, 00513 bool canDestroyInput=false) 00514 { 00515 FFTWGlobalConfiguration::Lock(); 00516 fftw_plan_with_nthreads(threads); 00517 // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE 00518 // because FFTW_ESTIMATE guarantee to not destroy the input 00519 unsigned roflags = flags; 00520 if( ! (flags & FFTW_ESTIMATE) ) 00521 { 00522 roflags = flags | FFTW_WISDOM_ONLY; 00523 } 00524 PlanType plan = fftw_plan_dft_r2c(rank,n,in,out,roflags); 00525 if( plan == NULL ) 00526 { 00527 // no wisdom available for that plan 00528 if( canDestroyInput ) 00529 { 00530 // just create the plan 00531 plan = fftw_plan_dft_r2c(rank,n,in,out,flags); 00532 } 00533 else 00534 { 00535 // lets create a plan with a fake input to generate the wisdom 00536 int total = 1; 00537 for( int i=0; i<rank; i++ ) 00538 { 00539 total *= n[i]; 00540 } 00541 PixelType * din = new PixelType[total]; 00542 fftw_plan_dft_r2c(rank,n,din,out,flags); 00543 delete [] din; 00544 // and then create the final plan - this time it shouldn't fail 00545 plan = fftw_plan_dft_r2c(rank,n,in,out,roflags); 00546 } 00547 FFTWGlobalConfiguration::SetNewWisdomAvailable(true); 00548 } 00549 FFTWGlobalConfiguration::Unlock(); 00550 assert( plan != NULL ); 00551 return plan; 00552 } 00553 00554 static PlanType Plan_dft_1d(int n, 00555 ComplexType *in, 00556 ComplexType *out, 00557 int sign, 00558 unsigned flags, 00559 int threads=1, 00560 bool canDestroyInput=false) 00561 { 00562 return Plan_dft(1, &n, in, out,sign , flags, threads, canDestroyInput); 00563 } 00564 00565 static PlanType Plan_dft_2d(int nx, 00566 int ny, 00567 ComplexType *in, 00568 ComplexType *out, 00569 int sign, 00570 unsigned flags, 00571 int threads=1, 00572 bool canDestroyInput=false) 00573 { 00574 int * sizes = new int[2]; 00575 sizes[0] = nx; 00576 sizes[1] = ny; 00577 PlanType plan = Plan_dft(2, sizes, in, out, sign, flags, threads, canDestroyInput); 00578 delete sizes; 00579 return plan; 00580 } 00581 00582 static PlanType Plan_dft_3d(int nx, 00583 int ny, 00584 int nz, 00585 ComplexType *in, 00586 ComplexType *out, 00587 int sign, 00588 unsigned flags, 00589 int threads=1, 00590 bool canDestroyInput=false) 00591 { 00592 int * sizes = new int[3]; 00593 sizes[0] = nx; 00594 sizes[1] = ny; 00595 sizes[2] = nz; 00596 PlanType plan = Plan_dft(3, sizes, in, out, sign, flags, threads, canDestroyInput); 00597 delete sizes; 00598 return plan; 00599 } 00600 00601 static PlanType Plan_dft(int rank, 00602 const int *n, 00603 ComplexType *in, 00604 ComplexType *out, 00605 int sign, 00606 unsigned flags, 00607 int threads=1, 00608 bool canDestroyInput=false) 00609 { 00610 FFTWGlobalConfiguration::Lock(); 00611 fftw_plan_with_nthreads(threads); 00612 // don't add FFTW_WISDOM_ONLY if the plan rigor is FFTW_ESTIMATE 00613 // because FFTW_ESTIMATE guarantee to not destroy the input 00614 unsigned roflags = flags; 00615 if( ! (flags & FFTW_ESTIMATE) ) 00616 { 00617 roflags = flags | FFTW_WISDOM_ONLY; 00618 } 00619 PlanType plan = fftw_plan_dft(rank,n,in,out,sign,roflags); 00620 if( plan == NULL ) 00621 { 00622 // no wisdom available for that plan 00623 if( canDestroyInput ) 00624 { 00625 // just create the plan 00626 plan = fftw_plan_dft(rank,n,in,out,sign,flags); 00627 } 00628 else 00629 { 00630 // lets create a plan with a fake input to generate the wisdom 00631 int total = 1; 00632 for( int i=0; i<rank; i++ ) 00633 { 00634 total *= n[i]; 00635 } 00636 ComplexType * din = new ComplexType[total]; 00637 fftw_plan_dft(rank,n,din,out,sign,flags); 00638 delete [] din; 00639 // and then create the final plan - this time it shouldn't fail 00640 plan = fftw_plan_dft(rank,n,in,out,sign,roflags); 00641 } 00642 FFTWGlobalConfiguration::SetNewWisdomAvailable(true); 00643 } 00644 FFTWGlobalConfiguration::Unlock(); 00645 assert( plan != NULL ); 00646 return plan; 00647 } 00648 00649 00650 static void Execute(PlanType p) 00651 { 00652 fftw_execute(p); 00653 } 00654 static void DestroyPlan(PlanType p) 00655 { 00656 fftw_destroy_plan(p); 00657 } 00658 }; 00659 00660 #endif 00661 } 00662 } 00663 #endif 00664
1.7.6.1