[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/multi_fft.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2009-2010 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_MULTI_FFT_HXX 00037 #define VIGRA_MULTI_FFT_HXX 00038 00039 #include "fftw3.hxx" 00040 #include "multi_array.hxx" 00041 #include "navigator.hxx" 00042 #include "copyimage.hxx" 00043 00044 namespace vigra { 00045 00046 /********************************************************/ 00047 /* */ 00048 /* Fourier Transform */ 00049 /* */ 00050 /********************************************************/ 00051 00052 /** \addtogroup FourierTransform 00053 */ 00054 //@{ 00055 00056 namespace detail { 00057 00058 template <unsigned int N, class T, class C> 00059 void moveDCToCenterImpl(MultiArrayView<N, T, C> a, unsigned int startDimension) 00060 { 00061 typedef typename MultiArrayView<N, T, C>::traverser Traverser; 00062 typedef MultiArrayNavigator<Traverser, N> Navigator; 00063 typedef typename Navigator::iterator Iterator; 00064 00065 for(unsigned int d = startDimension; d < N; ++d) 00066 { 00067 Navigator nav(a.traverser_begin(), a.shape(), d); 00068 00069 for( ; nav.hasMore(); nav++ ) 00070 { 00071 Iterator i = nav.begin(); 00072 int s = nav.end() - i; 00073 int s2 = s/2; 00074 00075 if(even(s)) 00076 { 00077 for(int k=0; k<s2; ++k) 00078 { 00079 std::swap(i[k], i[k+s2]); 00080 } 00081 } 00082 else 00083 { 00084 T v = i[0]; 00085 for(int k=0; k<s2; ++k) 00086 { 00087 i[k] = i[k+s2+1]; 00088 i[k+s2+1] = i[k+1]; 00089 } 00090 i[s2] = v; 00091 } 00092 } 00093 } 00094 } 00095 00096 template <unsigned int N, class T, class C> 00097 void moveDCToUpperLeftImpl(MultiArrayView<N, T, C> a, unsigned int startDimension) 00098 { 00099 typedef typename MultiArrayView<N, T, C>::traverser Traverser; 00100 typedef MultiArrayNavigator<Traverser, N> Navigator; 00101 typedef typename Navigator::iterator Iterator; 00102 00103 for(unsigned int d = startDimension; d < N; ++d) 00104 { 00105 Navigator nav(a.traverser_begin(), a.shape(), d); 00106 00107 for( ; nav.hasMore(); nav++ ) 00108 { 00109 Iterator i = nav.begin(); 00110 int s = nav.end() - i; 00111 int s2 = s/2; 00112 00113 if(even(s)) 00114 { 00115 for(int k=0; k<s2; ++k) 00116 { 00117 std::swap(i[k], i[k+s2]); 00118 } 00119 } 00120 else 00121 { 00122 T v = i[s2]; 00123 for(int k=s2; k>0; --k) 00124 { 00125 i[k] = i[k+s2]; 00126 i[k+s2] = i[k-1]; 00127 } 00128 i[0] = v; 00129 } 00130 } 00131 } 00132 } 00133 00134 } // namespace detail 00135 00136 /********************************************************/ 00137 /* */ 00138 /* moveDCToCenter */ 00139 /* */ 00140 /********************************************************/ 00141 00142 template <unsigned int N, class T, class C> 00143 inline void moveDCToCenter(MultiArrayView<N, T, C> a) 00144 { 00145 detail::moveDCToCenterImpl(a, 0); 00146 } 00147 00148 template <unsigned int N, class T, class C> 00149 inline void moveDCToUpperLeft(MultiArrayView<N, T, C> a) 00150 { 00151 detail::moveDCToUpperLeftImpl(a, 0); 00152 } 00153 00154 template <unsigned int N, class T, class C> 00155 inline void moveDCToHalfspaceCenter(MultiArrayView<N, T, C> a) 00156 { 00157 detail::moveDCToCenterImpl(a, 1); 00158 } 00159 00160 template <unsigned int N, class T, class C> 00161 inline void moveDCToHalfspaceUpperLeft(MultiArrayView<N, T, C> a) 00162 { 00163 detail::moveDCToUpperLeftImpl(a, 1); 00164 } 00165 00166 namespace detail 00167 { 00168 00169 inline fftw_plan 00170 fftwPlanCreate(unsigned int N, int* shape, 00171 FFTWComplex<double> * in, int* instrides, int instep, 00172 FFTWComplex<double> * out, int* outstrides, int outstep, 00173 int sign, unsigned int planner_flags) 00174 { 00175 return fftw_plan_many_dft(N, shape, 1, 00176 (fftw_complex *)in, instrides, instep, 0, 00177 (fftw_complex *)out, outstrides, outstep, 0, 00178 sign, planner_flags); 00179 } 00180 00181 inline fftw_plan 00182 fftwPlanCreate(unsigned int N, int* shape, 00183 double * in, int* instrides, int instep, 00184 FFTWComplex<double> * out, int* outstrides, int outstep, 00185 int /*sign is ignored*/, unsigned int planner_flags) 00186 { 00187 return fftw_plan_many_dft_r2c(N, shape, 1, 00188 in, instrides, instep, 0, 00189 (fftw_complex *)out, outstrides, outstep, 0, 00190 planner_flags); 00191 } 00192 00193 inline fftw_plan 00194 fftwPlanCreate(unsigned int N, int* shape, 00195 FFTWComplex<double> * in, int* instrides, int instep, 00196 double * out, int* outstrides, int outstep, 00197 int /*sign is ignored*/, unsigned int planner_flags) 00198 { 00199 return fftw_plan_many_dft_c2r(N, shape, 1, 00200 (fftw_complex *)in, instrides, instep, 0, 00201 out, outstrides, outstep, 0, 00202 planner_flags); 00203 } 00204 00205 inline fftwf_plan 00206 fftwPlanCreate(unsigned int N, int* shape, 00207 FFTWComplex<float> * in, int* instrides, int instep, 00208 FFTWComplex<float> * out, int* outstrides, int outstep, 00209 int sign, unsigned int planner_flags) 00210 { 00211 return fftwf_plan_many_dft(N, shape, 1, 00212 (fftwf_complex *)in, instrides, instep, 0, 00213 (fftwf_complex *)out, outstrides, outstep, 0, 00214 sign, planner_flags); 00215 } 00216 00217 inline fftwf_plan 00218 fftwPlanCreate(unsigned int N, int* shape, 00219 float * in, int* instrides, int instep, 00220 FFTWComplex<float> * out, int* outstrides, int outstep, 00221 int /*sign is ignored*/, unsigned int planner_flags) 00222 { 00223 return fftwf_plan_many_dft_r2c(N, shape, 1, 00224 in, instrides, instep, 0, 00225 (fftwf_complex *)out, outstrides, outstep, 0, 00226 planner_flags); 00227 } 00228 00229 inline fftwf_plan 00230 fftwPlanCreate(unsigned int N, int* shape, 00231 FFTWComplex<float> * in, int* instrides, int instep, 00232 float * out, int* outstrides, int outstep, 00233 int /*sign is ignored*/, unsigned int planner_flags) 00234 { 00235 return fftwf_plan_many_dft_c2r(N, shape, 1, 00236 (fftwf_complex *)in, instrides, instep, 0, 00237 out, outstrides, outstep, 0, 00238 planner_flags); 00239 } 00240 00241 inline fftwl_plan 00242 fftwPlanCreate(unsigned int N, int* shape, 00243 FFTWComplex<long double> * in, int* instrides, int instep, 00244 FFTWComplex<long double> * out, int* outstrides, int outstep, 00245 int sign, unsigned int planner_flags) 00246 { 00247 return fftwl_plan_many_dft(N, shape, 1, 00248 (fftwl_complex *)in, instrides, instep, 0, 00249 (fftwl_complex *)out, outstrides, outstep, 0, 00250 sign, planner_flags); 00251 } 00252 00253 inline fftwl_plan 00254 fftwPlanCreate(unsigned int N, int* shape, 00255 long double * in, int* instrides, int instep, 00256 FFTWComplex<long double> * out, int* outstrides, int outstep, 00257 int /*sign is ignored*/, unsigned int planner_flags) 00258 { 00259 return fftwl_plan_many_dft_r2c(N, shape, 1, 00260 in, instrides, instep, 0, 00261 (fftwl_complex *)out, outstrides, outstep, 0, 00262 planner_flags); 00263 } 00264 00265 inline fftwl_plan 00266 fftwPlanCreate(unsigned int N, int* shape, 00267 FFTWComplex<long double> * in, int* instrides, int instep, 00268 long double * out, int* outstrides, int outstep, 00269 int /*sign is ignored*/, unsigned int planner_flags) 00270 { 00271 return fftwl_plan_many_dft_c2r(N, shape, 1, 00272 (fftwl_complex *)in, instrides, instep, 0, 00273 out, outstrides, outstep, 0, 00274 planner_flags); 00275 } 00276 00277 inline void fftwPlanDestroy(fftw_plan plan) 00278 { 00279 if(plan != 0) 00280 fftw_destroy_plan(plan); 00281 } 00282 00283 inline void fftwPlanDestroy(fftwf_plan plan) 00284 { 00285 if(plan != 0) 00286 fftwf_destroy_plan(plan); 00287 } 00288 00289 inline void fftwPlanDestroy(fftwl_plan plan) 00290 { 00291 if(plan != 0) 00292 fftwl_destroy_plan(plan); 00293 } 00294 00295 inline void 00296 fftwPlanExecute(fftw_plan plan) 00297 { 00298 fftw_execute(plan); 00299 } 00300 00301 inline void 00302 fftwPlanExecute(fftwf_plan plan) 00303 { 00304 fftwf_execute(plan); 00305 } 00306 00307 inline void 00308 fftwPlanExecute(fftwl_plan plan) 00309 { 00310 fftwl_execute(plan); 00311 } 00312 00313 inline void 00314 fftwPlanExecute(fftw_plan plan, FFTWComplex<double> * in, FFTWComplex<double> * out) 00315 { 00316 fftw_execute_dft(plan, (fftw_complex *)in, (fftw_complex *)out); 00317 } 00318 00319 inline void 00320 fftwPlanExecute(fftw_plan plan, double * in, FFTWComplex<double> * out) 00321 { 00322 fftw_execute_dft_r2c(plan, in, (fftw_complex *)out); 00323 } 00324 00325 inline void 00326 fftwPlanExecute(fftw_plan plan, FFTWComplex<double> * in, double * out) 00327 { 00328 fftw_execute_dft_c2r(plan, (fftw_complex *)in, out); 00329 } 00330 00331 inline void 00332 fftwPlanExecute(fftwf_plan plan, FFTWComplex<float> * in, FFTWComplex<float> * out) 00333 { 00334 fftwf_execute_dft(plan, (fftwf_complex *)in, (fftwf_complex *)out); 00335 } 00336 00337 inline void 00338 fftwPlanExecute(fftwf_plan plan, float * in, FFTWComplex<float> * out) 00339 { 00340 fftwf_execute_dft_r2c(plan, in, (fftwf_complex *)out); 00341 } 00342 00343 inline void 00344 fftwPlanExecute(fftwf_plan plan, FFTWComplex<float> * in, float * out) 00345 { 00346 fftwf_execute_dft_c2r(plan, (fftwf_complex *)in, out); 00347 } 00348 00349 inline void 00350 fftwPlanExecute(fftwl_plan plan, FFTWComplex<long double> * in, FFTWComplex<long double> * out) 00351 { 00352 fftwl_execute_dft(plan, (fftwl_complex *)in, (fftwl_complex *)out); 00353 } 00354 00355 inline void 00356 fftwPlanExecute(fftwl_plan plan, long double * in, FFTWComplex<long double> * out) 00357 { 00358 fftwl_execute_dft_r2c(plan, in, (fftwl_complex *)out); 00359 } 00360 00361 inline void 00362 fftwPlanExecute(fftwl_plan plan, FFTWComplex<long double> * in, long double * out) 00363 { 00364 fftwl_execute_dft_c2r(plan, (fftwl_complex *)in, out); 00365 } 00366 00367 inline 00368 int fftwPaddingSize(int s) 00369 { 00370 // Image sizes where FFTW is fast. The list contains all 00371 // numbers less than 100000 whose prime decomposition is of the form 00372 // 2^a*3^b*5^c*7^d*11^e*13^f, where e+f is either 0 or 1, and the 00373 // other exponents are arbitrary 00374 static const int size = 1330; 00375 static int goodSizes[size] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 00376 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 35, 36, 39, 40, 42, 44, 45, 48, 00377 49, 50, 52, 54, 55, 56, 60, 63, 64, 65, 66, 70, 72, 75, 77, 78, 80, 81, 00378 84, 88, 90, 91, 96, 98, 99, 100, 104, 105, 108, 110, 112, 117, 120, 125, 00379 126, 128, 130, 132, 135, 140, 144, 147, 150, 154, 156, 160, 162, 165, 00380 168, 175, 176, 180, 182, 189, 192, 195, 196, 198, 200, 208, 210, 216, 00381 220, 224, 225, 231, 234, 240, 243, 245, 250, 252, 256, 260, 264, 270, 00382 273, 275, 280, 288, 294, 297, 300, 308, 312, 315, 320, 324, 325, 330, 00383 336, 343, 350, 351, 352, 360, 364, 375, 378, 384, 385, 390, 392, 396, 00384 400, 405, 416, 420, 432, 440, 441, 448, 450, 455, 462, 468, 480, 486, 00385 490, 495, 500, 504, 512, 520, 525, 528, 539, 540, 546, 550, 560, 567, 00386 576, 585, 588, 594, 600, 616, 624, 625, 630, 637, 640, 648, 650, 660, 00387 672, 675, 686, 693, 700, 702, 704, 720, 728, 729, 735, 750, 756, 768, 00388 770, 780, 784, 792, 800, 810, 819, 825, 832, 840, 864, 875, 880, 882, 00389 891, 896, 900, 910, 924, 936, 945, 960, 972, 975, 980, 990, 1000, 1008, 00390 1024, 1029, 1040, 1050, 1053, 1056, 1078, 1080, 1092, 1100, 1120, 1125, 00391 1134, 1152, 1155, 1170, 1176, 1188, 1200, 1215, 1225, 1232, 1248, 1250, 00392 1260, 1274, 1280, 1296, 1300, 1320, 1323, 1344, 1350, 1365, 1372, 1375, 00393 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470, 1485, 1500, 1512, 1536, 00394 1540, 1560, 1568, 1575, 1584, 1600, 1617, 1620, 1625, 1638, 1650, 1664, 00395 1680, 1701, 1715, 1728, 1750, 1755, 1760, 1764, 1782, 1792, 1800, 1820, 00396 1848, 1872, 1875, 1890, 1911, 1920, 1925, 1944, 1950, 1960, 1980, 2000, 00397 2016, 2025, 2048, 2058, 2079, 2080, 2100, 2106, 2112, 2156, 2160, 2184, 00398 2187, 2200, 2205, 2240, 2250, 2268, 2275, 2304, 2310, 2340, 2352, 2376, 00399 2400, 2401, 2430, 2450, 2457, 2464, 2475, 2496, 2500, 2520, 2548, 2560, 00400 2592, 2600, 2625, 2640, 2646, 2673, 2688, 2695, 2700, 2730, 2744, 2750, 00401 2772, 2800, 2808, 2816, 2835, 2880, 2912, 2916, 2925, 2940, 2970, 3000, 00402 3024, 3072, 3080, 3087, 3120, 3125, 3136, 3150, 3159, 3168, 3185, 3200, 00403 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3375, 3402, 3430, 3456, 3465, 00404 3500, 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3645, 3675, 3696, 3744, 00405 3750, 3773, 3780, 3822, 3840, 3850, 3888, 3900, 3920, 3960, 3969, 4000, 00406 4032, 4050, 4095, 4096, 4116, 4125, 4158, 4160, 4200, 4212, 4224, 4312, 00407 4320, 4368, 4374, 4375, 4400, 4410, 4455, 4459, 4480, 4500, 4536, 4550, 00408 4608, 4620, 4680, 4704, 4725, 4752, 4800, 4802, 4851, 4860, 4875, 4900, 00409 4914, 4928, 4950, 4992, 5000, 5040, 5096, 5103, 5120, 5145, 5184, 5200, 00410 5250, 5265, 5280, 5292, 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544, 00411 5600, 5616, 5625, 5632, 5670, 5733, 5760, 5775, 5824, 5832, 5850, 5880, 00412 5940, 6000, 6048, 6075, 6125, 6144, 6160, 6174, 6237, 6240, 6250, 6272, 00413 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6561, 6600, 6615, 00414 6656, 6720, 6750, 6804, 6825, 6860, 6875, 6912, 6930, 7000, 7020, 7040, 00415 7056, 7128, 7168, 7200, 7203, 7280, 7290, 7350, 7371, 7392, 7425, 7488, 00416 7500, 7546, 7560, 7644, 7680, 7700, 7776, 7800, 7840, 7875, 7920, 7938, 00417 8000, 8019, 8064, 8085, 8100, 8125, 8190, 8192, 8232, 8250, 8316, 8320, 00418 8400, 8424, 8448, 8505, 8575, 8624, 8640, 8736, 8748, 8750, 8775, 8800, 00419 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9261, 9360, 9375, 00420 9408, 9450, 9477, 9504, 9555, 9600, 9604, 9625, 9702, 9720, 9750, 9800, 00421 9828, 9856, 9900, 9984, 10000, 10080, 10125, 10192, 10206, 10240, 10290, 00422 10368, 10395, 10400, 10500, 10530, 10560, 10584, 10692, 10752, 10780, 00423 10800, 10920, 10935, 10976, 11000, 11025, 11088, 11200, 11232, 11250, 00424 11264, 11319, 11340, 11375, 11466, 11520, 11550, 11648, 11664, 11700, 00425 11760, 11880, 11907, 12000, 12005, 12096, 12150, 12250, 12285, 12288, 00426 12320, 12348, 12375, 12474, 12480, 12500, 12544, 12600, 12636, 12672, 00427 12740, 12800, 12936, 12960, 13000, 13104, 13122, 13125, 13200, 13230, 00428 13312, 13365, 13377, 13440, 13475, 13500, 13608, 13650, 13720, 13750, 00429 13824, 13860, 14000, 14040, 14080, 14112, 14175, 14256, 14336, 14400, 00430 14406, 14553, 14560, 14580, 14625, 14700, 14742, 14784, 14850, 14976, 00431 15000, 15092, 15120, 15288, 15309, 15360, 15400, 15435, 15552, 15600, 00432 15625, 15680, 15750, 15795, 15840, 15876, 15925, 16000, 16038, 16128, 00433 16170, 16200, 16250, 16380, 16384, 16464, 16500, 16632, 16640, 16800, 00434 16807, 16848, 16875, 16896, 17010, 17150, 17199, 17248, 17280, 17325, 00435 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000, 00436 18144, 18200, 18225, 18375, 18432, 18480, 18522, 18711, 18720, 18750, 00437 18816, 18865, 18900, 18954, 19008, 19110, 19200, 19208, 19250, 19404, 00438 19440, 19500, 19600, 19656, 19683, 19712, 19800, 19845, 19968, 20000, 00439 20160, 20250, 20384, 20412, 20475, 20480, 20580, 20625, 20736, 20790, 00440 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560, 21600, 21609, 00441 21840, 21870, 21875, 21952, 22000, 22050, 22113, 22176, 22275, 22295, 00442 22400, 22464, 22500, 22528, 22638, 22680, 22750, 22932, 23040, 23100, 00443 23296, 23328, 23400, 23520, 23625, 23760, 23814, 24000, 24010, 24057, 00444 24192, 24255, 24300, 24375, 24500, 24570, 24576, 24640, 24696, 24750, 00445 24948, 24960, 25000, 25088, 25200, 25272, 25344, 25480, 25515, 25600, 00446 25725, 25872, 25920, 26000, 26208, 26244, 26250, 26325, 26400, 26411, 00447 26460, 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440, 00448 27500, 27648, 27720, 27783, 28000, 28080, 28125, 28160, 28224, 28350, 00449 28431, 28512, 28665, 28672, 28800, 28812, 28875, 29106, 29120, 29160, 00450 29250, 29400, 29484, 29568, 29700, 29952, 30000, 30184, 30240, 30375, 00451 30576, 30618, 30625, 30720, 30800, 30870, 31104, 31185, 31200, 31213, 00452 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256, 00453 32340, 32400, 32500, 32760, 32768, 32805, 32928, 33000, 33075, 33264, 00454 33280, 33600, 33614, 33696, 33750, 33792, 33957, 34020, 34125, 34300, 00455 34375, 34398, 34496, 34560, 34650, 34944, 34992, 35000, 35100, 35200, 00456 35280, 35640, 35672, 35721, 35840, 36000, 36015, 36288, 36400, 36450, 00457 36750, 36855, 36864, 36960, 37044, 37125, 37422, 37440, 37500, 37632, 00458 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808, 38880, 00459 39000, 39200, 39312, 39366, 39375, 39424, 39600, 39690, 39936, 40000, 00460 40095, 40131, 40320, 40425, 40500, 40625, 40768, 40824, 40950, 40960, 00461 41160, 41250, 41472, 41580, 41600, 42000, 42120, 42240, 42336, 42525, 00462 42768, 42875, 43008, 43120, 43200, 43218, 43659, 43680, 43740, 43750, 00463 43875, 43904, 44000, 44100, 44226, 44352, 44550, 44590, 44800, 44928, 00464 45000, 45056, 45276, 45360, 45500, 45864, 45927, 46080, 46200, 46305, 00465 46592, 46656, 46800, 46875, 47040, 47250, 47385, 47520, 47628, 47775, 00466 48000, 48020, 48114, 48125, 48384, 48510, 48600, 48750, 49000, 49140, 00467 49152, 49280, 49392, 49500, 49896, 49920, 50000, 50176, 50400, 50421, 00468 50544, 50625, 50688, 50960, 51030, 51200, 51450, 51597, 51744, 51840, 00469 51975, 52000, 52416, 52488, 52500, 52650, 52800, 52822, 52920, 53248, 00470 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54675, 54880, 55000, 00471 55125, 55296, 55440, 55566, 56000, 56133, 56160, 56250, 56320, 56448, 00472 56595, 56700, 56862, 56875, 57024, 57330, 57344, 57600, 57624, 57750, 00473 58212, 58240, 58320, 58500, 58800, 58968, 59049, 59136, 59400, 59535, 00474 59904, 60000, 60025, 60368, 60480, 60750, 61152, 61236, 61250, 61425, 00475 61440, 61600, 61740, 61875, 62208, 62370, 62400, 62426, 62500, 62720, 00476 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512, 64680, 64800, 00477 64827, 65000, 65520, 65536, 65610, 65625, 65856, 66000, 66150, 66339, 00478 66528, 66560, 66825, 66885, 67200, 67228, 67375, 67392, 67500, 67584, 00479 67914, 68040, 68250, 68600, 68750, 68796, 68992, 69120, 69300, 69888, 00480 69984, 70000, 70200, 70400, 70560, 70875, 71280, 71344, 71442, 71680, 00481 72000, 72030, 72171, 72576, 72765, 72800, 72900, 73125, 73500, 73710, 00482 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264, 75460, 75600, 00483 75816, 76032, 76440, 76545, 76800, 76832, 77000, 77175, 77616, 77760, 00484 78000, 78125, 78400, 78624, 78732, 78750, 78848, 78975, 79200, 79233, 00485 79380, 79625, 79872, 80000, 80190, 80262, 80640, 80850, 81000, 81250, 00486 81536, 81648, 81900, 81920, 82320, 82500, 82944, 83160, 83200, 83349, 00487 84000, 84035, 84240, 84375, 84480, 84672, 85050, 85293, 85536, 85750, 00488 85995, 86016, 86240, 86400, 86436, 86625, 87318, 87360, 87480, 87500, 00489 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856, 00490 90000, 90112, 90552, 90720, 91000, 91125, 91728, 91854, 91875, 92160, 00491 92400, 92610, 93184, 93312, 93555, 93600, 93639, 93750, 94080, 94325, 00492 94500, 94770, 95040, 95256, 95550, 96000, 96040, 96228, 96250, 96768, 00493 97020, 97200, 97500, 98000, 98280, 98304, 98415, 98560, 98784, 99000, 00494 99225, 99792, 99840 }; 00495 00496 if(s <= 0 || s > goodSizes[size-1]) 00497 return s; 00498 return *std::upper_bound(goodSizes, goodSizes+size, s, std::less_equal<int>()); 00499 } 00500 00501 inline 00502 int fftwEvenPaddingSize(int s) 00503 { 00504 // Image sizes where FFTW is fast. The list contains all even 00505 // numbers less than 100000 whose prime decomposition is of the form 00506 // 2^a*3^b*5^c*7^d*11^e*13^f, where a >= 1, e+f is either 0 or 1, and the 00507 // other exponents are arbitrary 00508 static const int size = 1063; 00509 static int goodSizes[size] = { 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 00510 24, 26, 28, 30, 32, 36, 40, 42, 44, 48, 50, 52, 54, 56, 60, 64, 66, 70, 00511 72, 78, 80, 84, 88, 90, 96, 98, 100, 104, 108, 110, 112, 120, 126, 128, 00512 130, 132, 140, 144, 150, 154, 156, 160, 162, 168, 176, 180, 182, 192, 00513 196, 198, 200, 208, 210, 216, 220, 224, 234, 240, 250, 252, 256, 260, 00514 264, 270, 280, 288, 294, 300, 308, 312, 320, 324, 330, 336, 350, 352, 00515 360, 364, 378, 384, 390, 392, 396, 400, 416, 420, 432, 440, 448, 450, 00516 462, 468, 480, 486, 490, 500, 504, 512, 520, 528, 540, 546, 550, 560, 00517 576, 588, 594, 600, 616, 624, 630, 640, 648, 650, 660, 672, 686, 700, 00518 702, 704, 720, 728, 750, 756, 768, 770, 780, 784, 792, 800, 810, 832, 00519 840, 864, 880, 882, 896, 900, 910, 924, 936, 960, 972, 980, 990, 1000, 00520 1008, 1024, 1040, 1050, 1056, 1078, 1080, 1092, 1100, 1120, 1134, 1152, 00521 1170, 1176, 1188, 1200, 1232, 1248, 1250, 1260, 1274, 1280, 1296, 1300, 00522 1320, 1344, 1350, 1372, 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470, 00523 1500, 1512, 1536, 1540, 1560, 1568, 1584, 1600, 1620, 1638, 1650, 1664, 00524 1680, 1728, 1750, 1760, 1764, 1782, 1792, 1800, 1820, 1848, 1872, 1890, 00525 1920, 1944, 1950, 1960, 1980, 2000, 2016, 2048, 2058, 2080, 2100, 2106, 00526 2112, 2156, 2160, 2184, 2200, 2240, 2250, 2268, 2304, 2310, 2340, 2352, 00527 2376, 2400, 2430, 2450, 2464, 2496, 2500, 2520, 2548, 2560, 2592, 2600, 00528 2640, 2646, 2688, 2700, 2730, 2744, 2750, 2772, 2800, 2808, 2816, 2880, 00529 2912, 2916, 2940, 2970, 3000, 3024, 3072, 3080, 3120, 3136, 3150, 3168, 00530 3200, 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3402, 3430, 3456, 3500, 00531 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3696, 3744, 3750, 3780, 3822, 00532 3840, 3850, 3888, 3900, 3920, 3960, 4000, 4032, 4050, 4096, 4116, 4158, 00533 4160, 4200, 4212, 4224, 4312, 4320, 4368, 4374, 4400, 4410, 4480, 4500, 00534 4536, 4550, 4608, 4620, 4680, 4704, 4752, 4800, 4802, 4860, 4900, 4914, 00535 4928, 4950, 4992, 5000, 5040, 5096, 5120, 5184, 5200, 5250, 5280, 5292, 00536 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544, 5600, 5616, 5632, 5670, 00537 5760, 5824, 5832, 5850, 5880, 5940, 6000, 6048, 6144, 6160, 6174, 6240, 00538 6250, 6272, 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6600, 00539 6656, 6720, 6750, 6804, 6860, 6912, 6930, 7000, 7020, 7040, 7056, 7128, 00540 7168, 7200, 7280, 7290, 7350, 7392, 7488, 7500, 7546, 7560, 7644, 7680, 00541 7700, 7776, 7800, 7840, 7920, 7938, 8000, 8064, 8100, 8190, 8192, 8232, 00542 8250, 8316, 8320, 8400, 8424, 8448, 8624, 8640, 8736, 8748, 8750, 8800, 00543 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9360, 9408, 9450, 00544 9504, 9600, 9604, 9702, 9720, 9750, 9800, 9828, 9856, 9900, 9984, 10000, 00545 10080, 10192, 10206, 10240, 10290, 10368, 10400, 10500, 10530, 10560, 00546 10584, 10692, 10752, 10780, 10800, 10920, 10976, 11000, 11088, 11200, 00547 11232, 11250, 11264, 11340, 11466, 11520, 11550, 11648, 11664, 11700, 00548 11760, 11880, 12000, 12096, 12150, 12250, 12288, 12320, 12348, 12474, 00549 12480, 12500, 12544, 12600, 12636, 12672, 12740, 12800, 12936, 12960, 00550 13000, 13104, 13122, 13200, 13230, 13312, 13440, 13500, 13608, 13650, 00551 13720, 13750, 13824, 13860, 14000, 14040, 14080, 14112, 14256, 14336, 00552 14400, 14406, 14560, 14580, 14700, 14742, 14784, 14850, 14976, 15000, 00553 15092, 15120, 15288, 15360, 15400, 15552, 15600, 15680, 15750, 15840, 00554 15876, 16000, 16038, 16128, 16170, 16200, 16250, 16380, 16384, 16464, 00555 16500, 16632, 16640, 16800, 16848, 16896, 17010, 17150, 17248, 17280, 00556 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000, 00557 18144, 18200, 18432, 18480, 18522, 18720, 18750, 18816, 18900, 18954, 00558 19008, 19110, 19200, 19208, 19250, 19404, 19440, 19500, 19600, 19656, 00559 19712, 19800, 19968, 20000, 20160, 20250, 20384, 20412, 20480, 20580, 00560 20736, 20790, 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560, 00561 21600, 21840, 21870, 21952, 22000, 22050, 22176, 22400, 22464, 22500, 00562 22528, 22638, 22680, 22750, 22932, 23040, 23100, 23296, 23328, 23400, 00563 23520, 23760, 23814, 24000, 24010, 24192, 24300, 24500, 24570, 24576, 00564 24640, 24696, 24750, 24948, 24960, 25000, 25088, 25200, 25272, 25344, 00565 25480, 25600, 25872, 25920, 26000, 26208, 26244, 26250, 26400, 26460, 00566 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440, 27500, 00567 27648, 27720, 28000, 28080, 28160, 28224, 28350, 28512, 28672, 28800, 00568 28812, 29106, 29120, 29160, 29250, 29400, 29484, 29568, 29700, 29952, 00569 30000, 30184, 30240, 30576, 30618, 30720, 30800, 30870, 31104, 31200, 00570 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256, 00571 32340, 32400, 32500, 32760, 32768, 32928, 33000, 33264, 33280, 33600, 00572 33614, 33696, 33750, 33792, 34020, 34300, 34398, 34496, 34560, 34650, 00573 34944, 34992, 35000, 35100, 35200, 35280, 35640, 35672, 35840, 36000, 00574 36288, 36400, 36450, 36750, 36864, 36960, 37044, 37422, 37440, 37500, 00575 37632, 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808, 00576 38880, 39000, 39200, 39312, 39366, 39424, 39600, 39690, 39936, 40000, 00577 40320, 40500, 40768, 40824, 40950, 40960, 41160, 41250, 41472, 41580, 00578 41600, 42000, 42120, 42240, 42336, 42768, 43008, 43120, 43200, 43218, 00579 43680, 43740, 43750, 43904, 44000, 44100, 44226, 44352, 44550, 44590, 00580 44800, 44928, 45000, 45056, 45276, 45360, 45500, 45864, 46080, 46200, 00581 46592, 46656, 46800, 47040, 47250, 47520, 47628, 48000, 48020, 48114, 00582 48384, 48510, 48600, 48750, 49000, 49140, 49152, 49280, 49392, 49500, 00583 49896, 49920, 50000, 50176, 50400, 50544, 50688, 50960, 51030, 51200, 00584 51450, 51744, 51840, 52000, 52416, 52488, 52500, 52650, 52800, 52822, 00585 52920, 53248, 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54880, 00586 55000, 55296, 55440, 55566, 56000, 56160, 56250, 56320, 56448, 56700, 00587 56862, 57024, 57330, 57344, 57600, 57624, 57750, 58212, 58240, 58320, 00588 58500, 58800, 58968, 59136, 59400, 59904, 60000, 60368, 60480, 60750, 00589 61152, 61236, 61250, 61440, 61600, 61740, 62208, 62370, 62400, 62426, 00590 62500, 62720, 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512, 00591 64680, 64800, 65000, 65520, 65536, 65610, 65856, 66000, 66150, 66528, 00592 66560, 67200, 67228, 67392, 67500, 67584, 67914, 68040, 68250, 68600, 00593 68750, 68796, 68992, 69120, 69300, 69888, 69984, 70000, 70200, 70400, 00594 70560, 71280, 71344, 71442, 71680, 72000, 72030, 72576, 72800, 72900, 00595 73500, 73710, 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264, 00596 75460, 75600, 75816, 76032, 76440, 76800, 76832, 77000, 77616, 77760, 00597 78000, 78400, 78624, 78732, 78750, 78848, 79200, 79380, 79872, 80000, 00598 80190, 80262, 80640, 80850, 81000, 81250, 81536, 81648, 81900, 81920, 00599 82320, 82500, 82944, 83160, 83200, 84000, 84240, 84480, 84672, 85050, 00600 85536, 85750, 86016, 86240, 86400, 86436, 87318, 87360, 87480, 87500, 00601 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856, 00602 90000, 90112, 90552, 90720, 91000, 91728, 91854, 92160, 92400, 92610, 00603 93184, 93312, 93600, 93750, 94080, 94500, 94770, 95040, 95256, 95550, 00604 96000, 96040, 96228, 96250, 96768, 97020, 97200, 97500, 98000, 98280, 00605 98304, 98560, 98784, 99000, 99792, 99840 }; 00606 00607 if(s <= 0 || s > goodSizes[size-1]) 00608 return s; 00609 return *std::upper_bound(goodSizes, goodSizes+size, s, std::less_equal<int>()); 00610 } 00611 00612 template <int M> 00613 struct FFTEmbedKernel 00614 { 00615 template <unsigned int N, class Real, class C, class Shape> 00616 static void 00617 exec(MultiArrayView<N, Real, C> & out, Shape const & kernelShape, 00618 Shape & srcPoint, Shape & destPoint, bool copyIt) 00619 { 00620 for(srcPoint[M]=0; srcPoint[M]<kernelShape[M]; ++srcPoint[M]) 00621 { 00622 if(srcPoint[M] < (kernelShape[M] + 1) / 2) 00623 { 00624 destPoint[M] = srcPoint[M]; 00625 } 00626 else 00627 { 00628 destPoint[M] = srcPoint[M] + out.shape(M) - kernelShape[M]; 00629 copyIt = true; 00630 } 00631 FFTEmbedKernel<M-1>::exec(out, kernelShape, srcPoint, destPoint, copyIt); 00632 } 00633 } 00634 }; 00635 00636 template <> 00637 struct FFTEmbedKernel<0> 00638 { 00639 template <unsigned int N, class Real, class C, class Shape> 00640 static void 00641 exec(MultiArrayView<N, Real, C> & out, Shape const & kernelShape, 00642 Shape & srcPoint, Shape & destPoint, bool copyIt) 00643 { 00644 for(srcPoint[0]=0; srcPoint[0]<kernelShape[0]; ++srcPoint[0]) 00645 { 00646 if(srcPoint[0] < (kernelShape[0] + 1) / 2) 00647 { 00648 destPoint[0] = srcPoint[0]; 00649 } 00650 else 00651 { 00652 destPoint[0] = srcPoint[0] + out.shape(0) - kernelShape[0]; 00653 copyIt = true; 00654 } 00655 if(copyIt) 00656 { 00657 out[destPoint] = out[srcPoint]; 00658 out[srcPoint] = 0.0; 00659 } 00660 } 00661 } 00662 }; 00663 00664 template <unsigned int N, class Real, class C1, class C2> 00665 void 00666 fftEmbedKernel(MultiArrayView<N, Real, C1> kernel, 00667 MultiArrayView<N, Real, C2> out, 00668 Real norm = 1.0) 00669 { 00670 typedef typename MultiArrayShape<N>::type Shape; 00671 00672 MultiArrayView<N, Real, C2> kout = out.subarray(Shape(), kernel.shape()); 00673 00674 out.init(0.0); 00675 kout = kernel; 00676 if (norm != 1.0) 00677 kout *= norm; 00678 moveDCToUpperLeft(kout); 00679 00680 Shape srcPoint, destPoint; 00681 FFTEmbedKernel<(int)N-1>::exec(out, kernel.shape(), srcPoint, destPoint, false); 00682 } 00683 00684 template <unsigned int N, class Real, class C1, class C2> 00685 void 00686 fftEmbedArray(MultiArrayView<N, Real, C1> in, 00687 MultiArrayView<N, Real, C2> out) 00688 { 00689 typedef typename MultiArrayShape<N>::type Shape; 00690 00691 Shape diff = out.shape() - in.shape(), 00692 leftDiff = div(diff, MultiArrayIndex(2)), 00693 rightDiff = diff - leftDiff, 00694 right = in.shape() + leftDiff; 00695 00696 out.subarray(leftDiff, right) = in; 00697 00698 typedef typename MultiArrayView<N, Real, C2>::traverser Traverser; 00699 typedef MultiArrayNavigator<Traverser, N> Navigator; 00700 typedef typename Navigator::iterator Iterator; 00701 00702 for(unsigned int d = 0; d < N; ++d) 00703 { 00704 Navigator nav(out.traverser_begin(), out.shape(), d); 00705 00706 for( ; nav.hasMore(); nav++ ) 00707 { 00708 Iterator i = nav.begin(); 00709 for(int k=1; k<=leftDiff[d]; ++k) 00710 i[leftDiff[d] - k] = i[leftDiff[d] + k]; 00711 for(int k=0; k<rightDiff[d]; ++k) 00712 i[right[d] + k] = i[right[d] - k - 2]; 00713 } 00714 } 00715 } 00716 00717 } // namespace detail 00718 00719 template <class T, int N> 00720 TinyVector<T, N> 00721 fftwBestPaddedShape(TinyVector<T, N> shape) 00722 { 00723 for(unsigned int k=0; k<N; ++k) 00724 shape[k] = detail::fftwPaddingSize(shape[k]); 00725 return shape; 00726 } 00727 00728 template <class T, int N> 00729 TinyVector<T, N> 00730 fftwBestPaddedShapeR2C(TinyVector<T, N> shape) 00731 { 00732 shape[0] = detail::fftwEvenPaddingSize(shape[0]); 00733 for(unsigned int k=1; k<N; ++k) 00734 shape[k] = detail::fftwPaddingSize(shape[k]); 00735 return shape; 00736 } 00737 00738 /** \brief Find frequency domain shape for a R2C Fourier transform. 00739 00740 When a real valued array is transformed to the frequency domain, about half of the 00741 Fourier coefficients are redundant. The transform can be optimized as a <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C 00742 transform</a> that doesn't compute and store the redundant coefficients. This function 00743 computes the appropriate frequency domain shape for a given shape in the spatial domain. 00744 It simply replaces <tt>shape[0]</tt> with <tt>shape[0] / 2 + 1</tt>. 00745 00746 <b>\#include</b> <vigra/multi_fft.hxx> 00747 */ 00748 template <class T, int N> 00749 TinyVector<T, N> 00750 fftwCorrespondingShapeR2C(TinyVector<T, N> shape) 00751 { 00752 shape[0] = shape[0] / 2 + 1; 00753 return shape; 00754 } 00755 00756 template <class T, int N> 00757 TinyVector<T, N> 00758 fftwCorrespondingShapeC2R(TinyVector<T, N> shape, bool oddDimension0 = false) 00759 { 00760 shape[0] = oddDimension0 00761 ? (shape[0] - 1) * 2 + 1 00762 : (shape[0] - 1) * 2; 00763 return shape; 00764 } 00765 00766 /********************************************************/ 00767 /* */ 00768 /* FFTWPlan */ 00769 /* */ 00770 /********************************************************/ 00771 00772 /** C++ wrapper for FFTW plans. 00773 00774 The class encapsulates the calls to <tt>fftw_plan_dft_2d</tt>, <tt>fftw_execute</tt>, and 00775 <tt>fftw_destroy_plan</tt> (and their <tt>float</tt> and <tt>long double</tt> counterparts) 00776 in an easy-to-use interface. 00777 00778 Usually, you use this class only indirectly via \ref fourierTransform() 00779 and \ref fourierTransformInverse(). You only need this class if you want to have more control 00780 about FFTW's planning process (by providing non-default planning flags) and/or want to re-use 00781 plans for several transformations. 00782 00783 <b> Usage:</b> 00784 00785 <b>\#include</b> <vigra/multi_fft.hxx><br> 00786 Namespace: vigra 00787 00788 \code 00789 // compute complex Fourier transform of a real image 00790 MultiArray<2, double> src(Shape2(w, h)); 00791 MultiArray<2, FFTWComplex<double> > fourier(Shape2(w, h)); 00792 00793 // create an optimized plan by measuring the speed of several algorithm variants 00794 FFTWPlan<2, double> plan(src, fourier, FFTW_MEASURE); 00795 00796 plan.execute(src, fourier); 00797 \endcode 00798 */ 00799 template <unsigned int N, class Real = double> 00800 class FFTWPlan 00801 { 00802 typedef ArrayVector<int> Shape; 00803 typedef typename FFTWReal2Complex<Real>::plan_type PlanType; 00804 typedef typename FFTWComplex<Real>::complex_type Complex; 00805 00806 PlanType plan; 00807 Shape shape, instrides, outstrides; 00808 int sign; 00809 00810 public: 00811 /** \brief Create an empty plan. 00812 00813 The plan can be initialized later by one of the init() functions. 00814 */ 00815 FFTWPlan() 00816 : plan(0) 00817 {} 00818 00819 /** \brief Create a plan for a complex-to-complex transform. 00820 00821 \arg SIGN must be <tt>FFTW_FORWARD</tt> or <tt>FFTW_BACKWARD</tt> according to the 00822 desired transformation direction. 00823 \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 00824 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 00825 optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 00826 */ 00827 template <class C1, class C2> 00828 FFTWPlan(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00829 MultiArrayView<N, FFTWComplex<Real>, C2> out, 00830 int SIGN, unsigned int planner_flags = FFTW_ESTIMATE) 00831 : plan(0) 00832 { 00833 init(in, out, SIGN, planner_flags); 00834 } 00835 00836 /** \brief Create a plan for a real-to-complex transform. 00837 00838 This always refers to a forward transform. The shape of the output determines 00839 if a standard transform (when <tt>out.shape() == in.shape()</tt>) or an 00840 <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C 00841 transform</a> (when <tt>out.shape() == fftwCorrespondingShapeR2C(in.shape())</tt>) will be executed. 00842 00843 \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 00844 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 00845 optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 00846 */ 00847 template <class C1, class C2> 00848 FFTWPlan(MultiArrayView<N, Real, C1> in, 00849 MultiArrayView<N, FFTWComplex<Real>, C2> out, 00850 unsigned int planner_flags = FFTW_ESTIMATE) 00851 : plan(0) 00852 { 00853 init(in, out, planner_flags); 00854 } 00855 00856 /** \brief Create a plan for a complex-to-real transform. 00857 00858 This always refers to a inverse transform. The shape of the input determines 00859 if a standard transform (when <tt>in.shape() == out.shape()</tt>) or a 00860 <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">C2R 00861 transform</a> (when <tt>in.shape() == fftwCorrespondingShapeR2C(out.shape())</tt>) will be executed. 00862 00863 \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 00864 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 00865 optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 00866 */ 00867 template <class C1, class C2> 00868 FFTWPlan(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00869 MultiArrayView<N, Real, C2> out, 00870 unsigned int planner_flags = FFTW_ESTIMATE) 00871 : plan(0) 00872 { 00873 init(in, out, planner_flags); 00874 } 00875 00876 /** \brief Copy constructor. 00877 */ 00878 FFTWPlan(FFTWPlan const & other) 00879 : plan(other.plan), 00880 sign(other.sign) 00881 { 00882 FFTWPlan & o = const_cast<FFTWPlan &>(other); 00883 shape.swap(o.shape); 00884 instrides.swap(o.instrides); 00885 outstrides.swap(o.outstrides); 00886 o.plan = 0; // act like std::auto_ptr 00887 } 00888 00889 /** \brief Copy assigment. 00890 */ 00891 FFTWPlan & operator=(FFTWPlan const & other) 00892 { 00893 if(this != &other) 00894 { 00895 FFTWPlan & o = const_cast<FFTWPlan &>(other); 00896 plan = o.plan; 00897 shape.swap(o.shape); 00898 instrides.swap(o.instrides); 00899 outstrides.swap(o.outstrides); 00900 sign = o.sign; 00901 o.plan = 0; // act like std::auto_ptr 00902 } 00903 return *this; 00904 } 00905 00906 /** \brief Destructor. 00907 */ 00908 ~FFTWPlan() 00909 { 00910 detail::fftwPlanDestroy(plan); 00911 } 00912 00913 /** \brief Init a complex-to-complex transform. 00914 00915 See the constructor with the same signature for details. 00916 */ 00917 template <class C1, class C2> 00918 void init(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00919 MultiArrayView<N, FFTWComplex<Real>, C2> out, 00920 int SIGN, unsigned int planner_flags = FFTW_ESTIMATE) 00921 { 00922 vigra_precondition(in.strideOrdering() == out.strideOrdering(), 00923 "FFTWPlan.init(): input and output must have the same stride ordering."); 00924 00925 initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(), 00926 SIGN, planner_flags); 00927 } 00928 00929 /** \brief Init a real-to-complex transform. 00930 00931 See the constructor with the same signature for details. 00932 */ 00933 template <class C1, class C2> 00934 void init(MultiArrayView<N, Real, C1> in, 00935 MultiArrayView<N, FFTWComplex<Real>, C2> out, 00936 unsigned int planner_flags = FFTW_ESTIMATE) 00937 { 00938 vigra_precondition(in.strideOrdering() == out.strideOrdering(), 00939 "FFTWPlan.init(): input and output must have the same stride ordering."); 00940 00941 initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(), 00942 FFTW_FORWARD, planner_flags); 00943 } 00944 00945 /** \brief Init a complex-to-real transform. 00946 00947 See the constructor with the same signature for details. 00948 */ 00949 template <class C1, class C2> 00950 void init(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00951 MultiArrayView<N, Real, C2> out, 00952 unsigned int planner_flags = FFTW_ESTIMATE) 00953 { 00954 vigra_precondition(in.strideOrdering() == out.strideOrdering(), 00955 "FFTWPlan.init(): input and output must have the same stride ordering."); 00956 00957 initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(), 00958 FFTW_BACKWARD, planner_flags); 00959 } 00960 00961 /** \brief Execute a complex-to-complex transform. 00962 00963 The array shapes must be the same as in the corresponding init function 00964 or constructor. However, execute() can be called several times on 00965 the same plan, even with different arrays, as long as they have the appropriate 00966 shapes. 00967 */ 00968 template <class C1, class C2> 00969 void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00970 MultiArrayView<N, FFTWComplex<Real>, C2> out) const 00971 { 00972 executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending()); 00973 } 00974 00975 /** \brief Execute a real-to-complex transform. 00976 00977 The array shapes must be the same as in the corresponding init function 00978 or constructor. However, execute() can be called several times on 00979 the same plan, even with different arrays, as long as they have the appropriate 00980 shapes. 00981 */ 00982 template <class C1, class C2> 00983 void execute(MultiArrayView<N, Real, C1> in, 00984 MultiArrayView<N, FFTWComplex<Real>, C2> out) const 00985 { 00986 executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending()); 00987 } 00988 00989 /** \brief Execute a complex-to-real transform. 00990 00991 The array shapes must be the same as in the corresponding init function 00992 or constructor. However, execute() can be called several times on 00993 the same plan, even with different arrays, as long as they have the appropriate 00994 shapes. 00995 */ 00996 template <class C1, class C2> 00997 void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00998 MultiArrayView<N, Real, C2> out) const 00999 { 01000 executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending()); 01001 } 01002 01003 private: 01004 01005 template <class MI, class MO> 01006 void initImpl(MI ins, MO outs, int SIGN, unsigned int planner_flags); 01007 01008 template <class MI, class MO> 01009 void executeImpl(MI ins, MO outs) const; 01010 01011 void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> in, 01012 MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> out) const 01013 { 01014 vigra_precondition(in.shape() == out.shape(), 01015 "FFTWPlan.init(): input and output must have the same shape."); 01016 } 01017 01018 void checkShapes(MultiArrayView<N, Real, StridedArrayTag> ins, 01019 MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> outs) const 01020 { 01021 for(int k=0; k<(int)N-1; ++k) 01022 vigra_precondition(ins.shape(k) == outs.shape(k), 01023 "FFTWPlan.init(): input and output must have matching shapes."); 01024 vigra_precondition(ins.shape(N-1) / 2 + 1 == outs.shape(N-1), 01025 "FFTWPlan.init(): input and output must have matching shapes."); 01026 } 01027 01028 void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> ins, 01029 MultiArrayView<N, Real, StridedArrayTag> outs) const 01030 { 01031 for(int k=0; k<(int)N-1; ++k) 01032 vigra_precondition(ins.shape(k) == outs.shape(k), 01033 "FFTWPlan.init(): input and output must have matching shapes."); 01034 vigra_precondition(outs.shape(N-1) / 2 + 1 == ins.shape(N-1), 01035 "FFTWPlan.init(): input and output must have matching shapes."); 01036 } 01037 }; 01038 01039 template <unsigned int N, class Real> 01040 template <class MI, class MO> 01041 void 01042 FFTWPlan<N, Real>::initImpl(MI ins, MO outs, int SIGN, unsigned int planner_flags) 01043 { 01044 checkShapes(ins, outs); 01045 01046 typename MultiArrayShape<N>::type logicalShape(SIGN == FFTW_FORWARD 01047 ? ins.shape() 01048 : outs.shape()); 01049 01050 Shape newShape(logicalShape.begin(), logicalShape.end()), 01051 newIStrides(ins.stride().begin(), ins.stride().end()), 01052 newOStrides(outs.stride().begin(), outs.stride().end()), 01053 itotal(ins.shape().begin(), ins.shape().end()), 01054 ototal(outs.shape().begin(), outs.shape().end()); 01055 01056 for(unsigned int j=1; j<N; ++j) 01057 { 01058 itotal[j] = ins.stride(j-1) / ins.stride(j); 01059 ototal[j] = outs.stride(j-1) / outs.stride(j); 01060 } 01061 01062 PlanType newPlan = detail::fftwPlanCreate(N, newShape.begin(), 01063 ins.data(), itotal.begin(), ins.stride(N-1), 01064 outs.data(), ototal.begin(), outs.stride(N-1), 01065 SIGN, planner_flags); 01066 detail::fftwPlanDestroy(plan); 01067 plan = newPlan; 01068 shape.swap(newShape); 01069 instrides.swap(newIStrides); 01070 outstrides.swap(newOStrides); 01071 sign = SIGN; 01072 } 01073 01074 template <unsigned int N, class Real> 01075 template <class MI, class MO> 01076 void FFTWPlan<N, Real>::executeImpl(MI ins, MO outs) const 01077 { 01078 vigra_precondition(plan != 0, "FFTWPlan::execute(): plan is NULL."); 01079 01080 typename MultiArrayShape<N>::type lshape(sign == FFTW_FORWARD 01081 ? ins.shape() 01082 : outs.shape()); 01083 01084 vigra_precondition((lshape == TinyVectorView<int, N>(shape.data())), 01085 "FFTWPlan::execute(): shape mismatch between plan and data."); 01086 vigra_precondition((ins.stride() == TinyVectorView<int, N>(instrides.data())), 01087 "FFTWPlan::execute(): strides mismatch between plan and input data."); 01088 vigra_precondition((outs.stride() == TinyVectorView<int, N>(outstrides.data())), 01089 "FFTWPlan::execute(): strides mismatch between plan and output data."); 01090 01091 detail::fftwPlanExecute(plan, ins.data(), outs.data()); 01092 01093 typedef typename MO::value_type V; 01094 if(sign == FFTW_BACKWARD) 01095 outs *= V(1.0) / Real(outs.size()); 01096 } 01097 01098 /********************************************************/ 01099 /* */ 01100 /* FFTWConvolvePlan */ 01101 /* */ 01102 /********************************************************/ 01103 01104 /** C++ wrapper for a pair of FFTW plans used to perform FFT-based convolution. 01105 01106 The class encapsulates the calls to <tt>fftw_plan_dft_2d</tt>, <tt>fftw_execute</tt>, and 01107 <tt>fftw_destroy_plan</tt> (and their <tt>float</tt> and <tt>long double</tt> counterparts) 01108 in an easy-to-use interface. It always creates a pair of plans, one for the forward and one 01109 for the inverse transform required for convolution. 01110 01111 Usually, you use this class only indirectly via \ref convolveFFT() and its variants. 01112 You only need this class if you want to have more control about FFTW's planning process 01113 (by providing non-default planning flags) and/or want to re-use plans for several convolutions. 01114 01115 <b> Usage:</b> 01116 01117 <b>\#include</b> <vigra/multi_fft.hxx><br> 01118 Namespace: vigra 01119 01120 \code 01121 // convolve a real array with a real kernel 01122 MultiArray<2, double> src(Shape2(w, h)), dest(Shape2(w, h)); 01123 01124 MultiArray<2, double> spatial_kernel(Shape2(9, 9)); 01125 Gaussian<double> gauss(1.0); 01126 01127 for(int y=0; y<9; ++y) 01128 for(int x=0; x<9; ++x) 01129 spatial_kernel(x, y) = gauss(x-4.0)*gauss(y-4.0); 01130 01131 // create an optimized plan by measuring the speed of several algorithm variants 01132 FFTWConvolvePlan<2, double> plan(src, spatial_kernel, dest, FFTW_MEASURE); 01133 01134 plan.execute(src, spatial_kernel, dest); 01135 \endcode 01136 */ 01137 template <unsigned int N, class Real> 01138 class FFTWConvolvePlan 01139 { 01140 typedef FFTWComplex<Real> Complex; 01141 typedef MultiArrayView<N, Real, UnstridedArrayTag > RArray; 01142 typedef MultiArray<N, Complex, FFTWAllocator<Complex> > CArray; 01143 01144 FFTWPlan<N, Real> forward_plan, backward_plan; 01145 RArray realArray, realKernel; 01146 CArray fourierArray, fourierKernel; 01147 bool useFourierKernel; 01148 01149 public: 01150 01151 typedef typename MultiArrayShape<N>::type Shape; 01152 01153 /** \brief Create an empty plan. 01154 01155 The plan can be initialized later by one of the init() functions. 01156 */ 01157 FFTWConvolvePlan() 01158 : useFourierKernel(false) 01159 {} 01160 01161 /** \brief Create a plan to convolve a real array with a real kernel. 01162 01163 The kernel must be defined in the spatial domain. 01164 See \ref convolveFFT() for detailed information on required shapes and internal padding. 01165 01166 \arg planner_flags must be a combination of the 01167 <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 01168 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 01169 optimal algorithm settings or read them from pre-loaded 01170 <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 01171 */ 01172 template <class C1, class C2, class C3> 01173 FFTWConvolvePlan(MultiArrayView<N, Real, C1> in, 01174 MultiArrayView<N, Real, C2> kernel, 01175 MultiArrayView<N, Real, C3> out, 01176 unsigned int planner_flags = FFTW_ESTIMATE) 01177 : useFourierKernel(false) 01178 { 01179 init(in, kernel, out, planner_flags); 01180 } 01181 01182 /** \brief Create a plan to convolve a real array with a complex kernel. 01183 01184 The kernel must be defined in the Fourier domain, using the half-space format. 01185 See \ref convolveFFT() for detailed information on required shapes and internal padding. 01186 01187 \arg planner_flags must be a combination of the 01188 <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 01189 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 01190 optimal algorithm settings or read them from pre-loaded 01191 <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 01192 */ 01193 template <class C1, class C2, class C3> 01194 FFTWConvolvePlan(MultiArrayView<N, Real, C1> in, 01195 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01196 MultiArrayView<N, Real, C3> out, 01197 unsigned int planner_flags = FFTW_ESTIMATE) 01198 : useFourierKernel(true) 01199 { 01200 init(in, kernel, out, planner_flags); 01201 } 01202 01203 /** \brief Create a plan to convolve a complex array with a complex kernel. 01204 01205 See \ref convolveFFT() for detailed information on required shapes and internal padding. 01206 01207 \arg fourierDomainKernel determines if the kernel is defined in the spatial or 01208 Fourier domain. 01209 \arg planner_flags must be a combination of the 01210 <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 01211 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 01212 optimal algorithm settings or read them from pre-loaded 01213 <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 01214 */ 01215 template <class C1, class C2, class C3> 01216 FFTWConvolvePlan(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01217 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01218 MultiArrayView<N, FFTWComplex<Real>, C3> out, 01219 bool fourierDomainKernel, 01220 unsigned int planner_flags = FFTW_ESTIMATE) 01221 { 01222 init(in, kernel, out, fourierDomainKernel, planner_flags); 01223 } 01224 01225 01226 /** \brief Create a plan from just the shape information. 01227 01228 See \ref convolveFFT() for detailed information on required shapes and internal padding. 01229 01230 \arg fourierDomainKernel determines if the kernel is defined in the spatial or 01231 Fourier domain. 01232 \arg planner_flags must be a combination of the 01233 <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 01234 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 01235 optimal algorithm settings or read them from pre-loaded 01236 <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 01237 */ 01238 template <class C1, class C2, class C3> 01239 FFTWConvolvePlan(Shape inOut, Shape kernel, 01240 bool useFourierKernel = false, 01241 unsigned int planner_flags = FFTW_ESTIMATE) 01242 { 01243 if(useFourierKernel) 01244 init(inOut, kernel, planner_flags); 01245 else 01246 initFourierKernel(inOut, kernel, planner_flags); 01247 } 01248 01249 /** \brief Init a plan to convolve a real array with a real kernel. 01250 01251 See the constructor with the same signature for details. 01252 */ 01253 template <class C1, class C2, class C3> 01254 void init(MultiArrayView<N, Real, C1> in, 01255 MultiArrayView<N, Real, C2> kernel, 01256 MultiArrayView<N, Real, C3> out, 01257 unsigned int planner_flags = FFTW_ESTIMATE) 01258 { 01259 vigra_precondition(in.shape() == out.shape(), 01260 "FFTWConvolvePlan::init(): input and output must have the same shape."); 01261 init(in.shape(), kernel.shape(), planner_flags); 01262 } 01263 01264 /** \brief Init a plan to convolve a real array with a complex kernel. 01265 01266 See the constructor with the same signature for details. 01267 */ 01268 template <class C1, class C2, class C3> 01269 void init(MultiArrayView<N, Real, C1> in, 01270 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01271 MultiArrayView<N, Real, C3> out, 01272 unsigned int planner_flags = FFTW_ESTIMATE) 01273 { 01274 vigra_precondition(in.shape() == out.shape(), 01275 "FFTWConvolvePlan::init(): input and output must have the same shape."); 01276 initFourierKernel(in.shape(), kernel.shape(), planner_flags); 01277 } 01278 01279 /** \brief Init a plan to convolve a complex array with a complex kernel. 01280 01281 See the constructor with the same signature for details. 01282 */ 01283 template <class C1, class C2, class C3> 01284 void init(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01285 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01286 MultiArrayView<N, FFTWComplex<Real>, C3> out, 01287 bool fourierDomainKernel, 01288 unsigned int planner_flags = FFTW_ESTIMATE) 01289 { 01290 vigra_precondition(in.shape() == out.shape(), 01291 "FFTWConvolvePlan::init(): input and output must have the same shape."); 01292 useFourierKernel = fourierDomainKernel; 01293 initComplex(in.shape(), kernel.shape(), planner_flags); 01294 } 01295 01296 /** \brief Init a plan to convolve a real array with a sequence of kernels. 01297 01298 The kernels can be either real or complex. The sequences \a kernels and \a outs 01299 must have the same length. See the corresponding constructors 01300 for single kernels for details. 01301 */ 01302 template <class C1, class KernelIterator, class OutIterator> 01303 void initMany(MultiArrayView<N, Real, C1> in, 01304 KernelIterator kernels, KernelIterator kernelsEnd, 01305 OutIterator outs, unsigned int planner_flags = FFTW_ESTIMATE) 01306 { 01307 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray; 01308 typedef typename KernelArray::value_type KernelValue; 01309 typedef typename std::iterator_traits<OutIterator>::value_type OutArray; 01310 typedef typename OutArray::value_type OutValue; 01311 01312 bool realKernel = IsSameType<KernelValue, Real>::value; 01313 bool fourierKernel = IsSameType<KernelValue, Complex>::value; 01314 01315 vigra_precondition(realKernel || fourierKernel, 01316 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type."); 01317 vigra_precondition((IsSameType<OutValue, Real>::value), 01318 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type."); 01319 01320 if(realKernel) 01321 { 01322 initMany(in.shape(), checkShapes(in.shape(), kernels, kernelsEnd, outs), 01323 planner_flags); 01324 } 01325 else 01326 { 01327 initFourierKernelMany(in.shape(), 01328 checkShapesFourier(in.shape(), kernels, kernelsEnd, outs), 01329 planner_flags); 01330 } 01331 } 01332 01333 /** \brief Init a plan to convolve a complex array with a sequence of kernels. 01334 01335 The kernels must be complex as well. The sequences \a kernels and \a outs 01336 must have the same length. See the corresponding constructors 01337 for single kernels for details. 01338 */ 01339 template <class C1, class KernelIterator, class OutIterator> 01340 void initMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01341 KernelIterator kernels, KernelIterator kernelsEnd, 01342 OutIterator outs, 01343 bool fourierDomainKernels, 01344 unsigned int planner_flags = FFTW_ESTIMATE) 01345 { 01346 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray; 01347 typedef typename KernelArray::value_type KernelValue; 01348 typedef typename std::iterator_traits<OutIterator>::value_type OutArray; 01349 typedef typename OutArray::value_type OutValue; 01350 01351 vigra_precondition((IsSameType<KernelValue, Complex>::value), 01352 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type."); 01353 vigra_precondition((IsSameType<OutValue, Complex>::value), 01354 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type."); 01355 01356 useFourierKernel = fourierDomainKernels; 01357 01358 Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs); 01359 01360 CArray newFourierArray(paddedShape), newFourierKernel(paddedShape); 01361 01362 FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags); 01363 FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags); 01364 01365 forward_plan = fplan; 01366 backward_plan = bplan; 01367 fourierArray.swap(newFourierArray); 01368 fourierKernel.swap(newFourierKernel); 01369 } 01370 01371 void init(Shape inOut, Shape kernel, 01372 unsigned int planner_flags = FFTW_ESTIMATE); 01373 01374 void initFourierKernel(Shape inOut, Shape kernel, 01375 unsigned int planner_flags = FFTW_ESTIMATE); 01376 01377 void initComplex(Shape inOut, Shape kernel, 01378 unsigned int planner_flags = FFTW_ESTIMATE); 01379 01380 void initMany(Shape inOut, Shape maxKernel, 01381 unsigned int planner_flags = FFTW_ESTIMATE) 01382 { 01383 init(inOut, maxKernel, planner_flags); 01384 } 01385 01386 void initFourierKernelMany(Shape inOut, Shape kernels, 01387 unsigned int planner_flags = FFTW_ESTIMATE) 01388 { 01389 initFourierKernel(inOut, kernels, planner_flags); 01390 } 01391 01392 /** \brief Execute a plan to convolve a real array with a real kernel. 01393 01394 The array shapes must be the same as in the corresponding init function 01395 or constructor. However, execute() can be called several times on 01396 the same plan, even with different arrays, as long as they have the appropriate 01397 shapes. 01398 */ 01399 template <class C1, class C2, class C3> 01400 void execute(MultiArrayView<N, Real, C1> in, 01401 MultiArrayView<N, Real, C2> kernel, 01402 MultiArrayView<N, Real, C3> out); 01403 01404 /** \brief Execute a plan to convolve a real array with a complex kernel. 01405 01406 The array shapes must be the same as in the corresponding init function 01407 or constructor. However, execute() can be called several times on 01408 the same plan, even with different arrays, as long as they have the appropriate 01409 shapes. 01410 */ 01411 template <class C1, class C2, class C3> 01412 void execute(MultiArrayView<N, Real, C1> in, 01413 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01414 MultiArrayView<N, Real, C3> out); 01415 01416 /** \brief Execute a plan to convolve a complex array with a complex kernel. 01417 01418 The array shapes must be the same as in the corresponding init function 01419 or constructor. However, execute() can be called several times on 01420 the same plan, even with different arrays, as long as they have the appropriate 01421 shapes. 01422 */ 01423 template <class C1, class C2, class C3> 01424 void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01425 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01426 MultiArrayView<N, FFTWComplex<Real>, C3> out); 01427 01428 01429 /** \brief Execute a plan to convolve a complex array with a sequence of kernels. 01430 01431 The array shapes must be the same as in the corresponding init function 01432 or constructor. However, executeMany() can be called several times on 01433 the same plan, even with different arrays, as long as they have the appropriate 01434 shapes. 01435 */ 01436 template <class C1, class KernelIterator, class OutIterator> 01437 void executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01438 KernelIterator kernels, KernelIterator kernelsEnd, 01439 OutIterator outs); 01440 01441 /** \brief Execute a plan to convolve a real array with a sequence of kernels. 01442 01443 The array shapes must be the same as in the corresponding init function 01444 or constructor. However, executeMany() can be called several times on 01445 the same plan, even with different arrays, as long as they have the appropriate 01446 shapes. 01447 */ 01448 template <class C1, class KernelIterator, class OutIterator> 01449 void executeMany(MultiArrayView<N, Real, C1> in, 01450 KernelIterator kernels, KernelIterator kernelsEnd, 01451 OutIterator outs) 01452 { 01453 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray; 01454 typedef typename KernelArray::value_type KernelValue; 01455 typedef typename IsSameType<KernelValue, Complex>::type UseFourierKernel; 01456 typedef typename std::iterator_traits<OutIterator>::value_type OutArray; 01457 typedef typename OutArray::value_type OutValue; 01458 01459 bool realKernel = IsSameType<KernelValue, Real>::value; 01460 bool fourierKernel = IsSameType<KernelValue, Complex>::value; 01461 01462 vigra_precondition(realKernel || fourierKernel, 01463 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type."); 01464 vigra_precondition((IsSameType<OutValue, Real>::value), 01465 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type."); 01466 01467 executeManyImpl(in, kernels, kernelsEnd, outs, UseFourierKernel()); 01468 } 01469 01470 private: 01471 01472 template <class KernelIterator, class OutIterator> 01473 Shape checkShapes(Shape in, 01474 KernelIterator kernels, KernelIterator kernelsEnd, 01475 OutIterator outs); 01476 01477 template <class KernelIterator, class OutIterator> 01478 Shape checkShapesFourier(Shape in, 01479 KernelIterator kernels, KernelIterator kernelsEnd, 01480 OutIterator outs); 01481 01482 template <class KernelIterator, class OutIterator> 01483 Shape checkShapesComplex(Shape in, 01484 KernelIterator kernels, KernelIterator kernelsEnd, 01485 OutIterator outs); 01486 01487 template <class C1, class KernelIterator, class OutIterator> 01488 void 01489 executeManyImpl(MultiArrayView<N, Real, C1> in, 01490 KernelIterator kernels, KernelIterator kernelsEnd, 01491 OutIterator outs, VigraFalseType /* useFourierKernel*/); 01492 01493 template <class C1, class KernelIterator, class OutIterator> 01494 void 01495 executeManyImpl(MultiArrayView<N, Real, C1> in, 01496 KernelIterator kernels, KernelIterator kernelsEnd, 01497 OutIterator outs, VigraTrueType /* useFourierKernel*/); 01498 01499 }; 01500 01501 template <unsigned int N, class Real> 01502 void 01503 FFTWConvolvePlan<N, Real>::init(Shape in, Shape kernel, 01504 unsigned int planner_flags) 01505 { 01506 Shape paddedShape = fftwBestPaddedShapeR2C(in + kernel - Shape(1)), 01507 complexShape = fftwCorrespondingShapeR2C(paddedShape); 01508 01509 CArray newFourierArray(complexShape), newFourierKernel(complexShape); 01510 01511 Shape realStrides = 2*newFourierArray.stride(); 01512 realStrides[0] = 1; 01513 RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data()); 01514 RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data()); 01515 01516 FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags); 01517 FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags); 01518 01519 forward_plan = fplan; 01520 backward_plan = bplan; 01521 realArray = newRealArray; 01522 realKernel = newRealKernel; 01523 fourierArray.swap(newFourierArray); 01524 fourierKernel.swap(newFourierKernel); 01525 useFourierKernel = false; 01526 } 01527 01528 template <unsigned int N, class Real> 01529 void 01530 FFTWConvolvePlan<N, Real>::initFourierKernel(Shape in, Shape kernel, 01531 unsigned int planner_flags) 01532 { 01533 Shape complexShape = kernel, 01534 paddedShape = fftwCorrespondingShapeC2R(complexShape); 01535 01536 for(unsigned int k=0; k<N; ++k) 01537 vigra_precondition(in[k] <= paddedShape[k], 01538 "FFTWConvolvePlan::init(): kernel too small for given input."); 01539 01540 CArray newFourierArray(complexShape), newFourierKernel(complexShape); 01541 01542 Shape realStrides = 2*newFourierArray.stride(); 01543 realStrides[0] = 1; 01544 RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data()); 01545 RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data()); 01546 01547 FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags); 01548 FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags); 01549 01550 forward_plan = fplan; 01551 backward_plan = bplan; 01552 realArray = newRealArray; 01553 realKernel = newRealKernel; 01554 fourierArray.swap(newFourierArray); 01555 fourierKernel.swap(newFourierKernel); 01556 useFourierKernel = true; 01557 } 01558 01559 template <unsigned int N, class Real> 01560 void 01561 FFTWConvolvePlan<N, Real>::initComplex(Shape in, Shape kernel, 01562 unsigned int planner_flags) 01563 { 01564 Shape paddedShape; 01565 01566 if(useFourierKernel) 01567 { 01568 for(unsigned int k=0; k<N; ++k) 01569 vigra_precondition(in[k] <= kernel[k], 01570 "FFTWConvolvePlan::init(): kernel too small for given input."); 01571 01572 paddedShape = kernel; 01573 } 01574 else 01575 { 01576 paddedShape = fftwBestPaddedShape(in + kernel - Shape(1)); 01577 } 01578 01579 CArray newFourierArray(paddedShape), newFourierKernel(paddedShape); 01580 01581 FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags); 01582 FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags); 01583 01584 forward_plan = fplan; 01585 backward_plan = bplan; 01586 fourierArray.swap(newFourierArray); 01587 fourierKernel.swap(newFourierKernel); 01588 } 01589 01590 #ifndef DOXYGEN // doxygen documents these functions as free functions 01591 01592 template <unsigned int N, class Real> 01593 template <class C1, class C2, class C3> 01594 void 01595 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, Real, C1> in, 01596 MultiArrayView<N, Real, C2> kernel, 01597 MultiArrayView<N, Real, C3> out) 01598 { 01599 vigra_precondition(!useFourierKernel, 01600 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel."); 01601 01602 vigra_precondition(in.shape() == out.shape(), 01603 "FFTWConvolvePlan::execute(): input and output must have the same shape."); 01604 01605 Shape paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernel.shape() - Shape(1)), 01606 diff = paddedShape - in.shape(), 01607 left = div(diff, MultiArrayIndex(2)), 01608 right = in.shape() + left; 01609 01610 vigra_precondition(paddedShape == realArray.shape(), 01611 "FFTWConvolvePlan::execute(): shape mismatch between input and plan."); 01612 01613 detail::fftEmbedArray(in, realArray); 01614 forward_plan.execute(realArray, fourierArray); 01615 01616 detail::fftEmbedKernel(kernel, realKernel); 01617 forward_plan.execute(realKernel, fourierKernel); 01618 01619 fourierArray *= fourierKernel; 01620 01621 backward_plan.execute(fourierArray, realArray); 01622 01623 out = realArray.subarray(left, right); 01624 } 01625 01626 template <unsigned int N, class Real> 01627 template <class C1, class C2, class C3> 01628 void 01629 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, Real, C1> in, 01630 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01631 MultiArrayView<N, Real, C3> out) 01632 { 01633 vigra_precondition(useFourierKernel, 01634 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel."); 01635 01636 vigra_precondition(in.shape() == out.shape(), 01637 "FFTWConvolvePlan::execute(): input and output must have the same shape."); 01638 01639 vigra_precondition(kernel.shape() == fourierArray.shape(), 01640 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan."); 01641 01642 Shape paddedShape = fftwCorrespondingShapeC2R(kernel.shape(), odd(in.shape(0))), 01643 diff = paddedShape - in.shape(), 01644 left = div(diff, MultiArrayIndex(2)), 01645 right = in.shape() + left; 01646 01647 vigra_precondition(paddedShape == realArray.shape(), 01648 "FFTWConvolvePlan::execute(): shape mismatch between input and plan."); 01649 01650 detail::fftEmbedArray(in, realArray); 01651 forward_plan.execute(realArray, fourierArray); 01652 01653 fourierKernel = kernel; 01654 moveDCToHalfspaceUpperLeft(fourierKernel); 01655 01656 fourierArray *= fourierKernel; 01657 01658 backward_plan.execute(fourierArray, realArray); 01659 01660 out = realArray.subarray(left, right); 01661 } 01662 01663 template <unsigned int N, class Real> 01664 template <class C1, class C2, class C3> 01665 void 01666 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01667 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01668 MultiArrayView<N, FFTWComplex<Real>, C3> out) 01669 { 01670 vigra_precondition(in.shape() == out.shape(), 01671 "FFTWConvolvePlan::execute(): input and output must have the same shape."); 01672 01673 Shape paddedShape = fourierArray.shape(), 01674 diff = paddedShape - in.shape(), 01675 left = div(diff, MultiArrayIndex(2)), 01676 right = in.shape() + left; 01677 01678 if(useFourierKernel) 01679 { 01680 vigra_precondition(kernel.shape() == fourierArray.shape(), 01681 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan."); 01682 01683 fourierKernel = kernel; 01684 moveDCToUpperLeft(fourierKernel); 01685 } 01686 else 01687 { 01688 detail::fftEmbedKernel(kernel, fourierKernel); 01689 forward_plan.execute(fourierKernel, fourierKernel); 01690 } 01691 01692 detail::fftEmbedArray(in, fourierArray); 01693 forward_plan.execute(fourierArray, fourierArray); 01694 01695 fourierArray *= fourierKernel; 01696 01697 backward_plan.execute(fourierArray, fourierArray); 01698 01699 out = fourierArray.subarray(left, right); 01700 } 01701 01702 template <unsigned int N, class Real> 01703 template <class C1, class KernelIterator, class OutIterator> 01704 void 01705 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in, 01706 KernelIterator kernels, KernelIterator kernelsEnd, 01707 OutIterator outs, VigraFalseType /*useFourierKernel*/) 01708 { 01709 vigra_precondition(!useFourierKernel, 01710 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel."); 01711 01712 Shape kernelMax = checkShapes(in.shape(), kernels, kernelsEnd, outs), 01713 paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernelMax - Shape(1)), 01714 diff = paddedShape - in.shape(), 01715 left = div(diff, MultiArrayIndex(2)), 01716 right = in.shape() + left; 01717 01718 vigra_precondition(paddedShape == realArray.shape(), 01719 "FFTWConvolvePlan::executeMany(): shape mismatch between input and plan."); 01720 01721 detail::fftEmbedArray(in, realArray); 01722 forward_plan.execute(realArray, fourierArray); 01723 01724 for(; kernels != kernelsEnd; ++kernels, ++outs) 01725 { 01726 detail::fftEmbedKernel(*kernels, realKernel); 01727 forward_plan.execute(realKernel, fourierKernel); 01728 01729 fourierKernel *= fourierArray; 01730 01731 backward_plan.execute(fourierKernel, realKernel); 01732 01733 *outs = realKernel.subarray(left, right); 01734 } 01735 } 01736 01737 template <unsigned int N, class Real> 01738 template <class C1, class KernelIterator, class OutIterator> 01739 void 01740 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in, 01741 KernelIterator kernels, KernelIterator kernelsEnd, 01742 OutIterator outs, VigraTrueType /*useFourierKernel*/) 01743 { 01744 vigra_precondition(useFourierKernel, 01745 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel."); 01746 01747 Shape complexShape = checkShapesFourier(in.shape(), kernels, kernelsEnd, outs), 01748 paddedShape = fftwCorrespondingShapeC2R(complexShape, odd(in.shape(0))), 01749 diff = paddedShape - in.shape(), 01750 left = div(diff, MultiArrayIndex(2)), 01751 right = in.shape() + left; 01752 01753 vigra_precondition(complexShape == fourierArray.shape(), 01754 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between kernels and plan."); 01755 01756 vigra_precondition(paddedShape == realArray.shape(), 01757 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between input and plan."); 01758 01759 detail::fftEmbedArray(in, realArray); 01760 forward_plan.execute(realArray, fourierArray); 01761 01762 for(; kernels != kernelsEnd; ++kernels, ++outs) 01763 { 01764 fourierKernel = *kernels; 01765 moveDCToHalfspaceUpperLeft(fourierKernel); 01766 fourierKernel *= fourierArray; 01767 01768 backward_plan.execute(fourierKernel, realKernel); 01769 01770 *outs = realKernel.subarray(left, right); 01771 } 01772 } 01773 01774 template <unsigned int N, class Real> 01775 template <class C1, class KernelIterator, class OutIterator> 01776 void 01777 FFTWConvolvePlan<N, Real>::executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01778 KernelIterator kernels, KernelIterator kernelsEnd, 01779 OutIterator outs) 01780 { 01781 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray; 01782 typedef typename KernelArray::value_type KernelValue; 01783 typedef typename std::iterator_traits<OutIterator>::value_type OutArray; 01784 typedef typename OutArray::value_type OutValue; 01785 01786 vigra_precondition((IsSameType<KernelValue, Complex>::value), 01787 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type."); 01788 vigra_precondition((IsSameType<OutValue, Complex>::value), 01789 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type."); 01790 01791 Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs), 01792 diff = paddedShape - in.shape(), 01793 left = div(diff, MultiArrayIndex(2)), 01794 right = in.shape() + left; 01795 01796 detail::fftEmbedArray(in, fourierArray); 01797 forward_plan.execute(fourierArray, fourierArray); 01798 01799 for(; kernels != kernelsEnd; ++kernels, ++outs) 01800 { 01801 if(useFourierKernel) 01802 { 01803 fourierKernel = *kernels; 01804 moveDCToUpperLeft(fourierKernel); 01805 } 01806 else 01807 { 01808 detail::fftEmbedKernel(*kernels, fourierKernel); 01809 forward_plan.execute(fourierKernel, fourierKernel); 01810 } 01811 01812 fourierKernel *= fourierArray; 01813 01814 backward_plan.execute(fourierKernel, fourierKernel); 01815 01816 *outs = fourierKernel.subarray(left, right); 01817 } 01818 } 01819 01820 #endif // DOXYGEN 01821 01822 template <unsigned int N, class Real> 01823 template <class KernelIterator, class OutIterator> 01824 typename FFTWConvolvePlan<N, Real>::Shape 01825 FFTWConvolvePlan<N, Real>::checkShapes(Shape in, 01826 KernelIterator kernels, KernelIterator kernelsEnd, 01827 OutIterator outs) 01828 { 01829 vigra_precondition(kernels != kernelsEnd, 01830 "FFTWConvolvePlan::checkShapes(): empty kernel sequence."); 01831 01832 Shape kernelMax; 01833 for(; kernels != kernelsEnd; ++kernels, ++outs) 01834 { 01835 vigra_precondition(in == outs->shape(), 01836 "FFTWConvolvePlan::checkShapes(): shape mismatch between input and (one) output."); 01837 kernelMax = max(kernelMax, kernels->shape()); 01838 } 01839 vigra_precondition(prod(kernelMax) > 0, 01840 "FFTWConvolvePlan::checkShapes(): all kernels have size 0."); 01841 return kernelMax; 01842 } 01843 01844 template <unsigned int N, class Real> 01845 template <class KernelIterator, class OutIterator> 01846 typename FFTWConvolvePlan<N, Real>::Shape 01847 FFTWConvolvePlan<N, Real>::checkShapesFourier(Shape in, 01848 KernelIterator kernels, KernelIterator kernelsEnd, 01849 OutIterator outs) 01850 { 01851 vigra_precondition(kernels != kernelsEnd, 01852 "FFTWConvolvePlan::checkShapesFourier(): empty kernel sequence."); 01853 01854 Shape complexShape = kernels->shape(), 01855 paddedShape = fftwCorrespondingShapeC2R(complexShape); 01856 01857 for(unsigned int k=0; k<N; ++k) 01858 vigra_precondition(in[k] <= paddedShape[k], 01859 "FFTWConvolvePlan::checkShapesFourier(): kernels too small for given input."); 01860 01861 for(; kernels != kernelsEnd; ++kernels, ++outs) 01862 { 01863 vigra_precondition(in == outs->shape(), 01864 "FFTWConvolvePlan::checkShapesFourier(): shape mismatch between input and (one) output."); 01865 vigra_precondition(complexShape == kernels->shape(), 01866 "FFTWConvolvePlan::checkShapesFourier(): all kernels must have the same size."); 01867 } 01868 return complexShape; 01869 } 01870 01871 template <unsigned int N, class Real> 01872 template <class KernelIterator, class OutIterator> 01873 typename FFTWConvolvePlan<N, Real>::Shape 01874 FFTWConvolvePlan<N, Real>::checkShapesComplex(Shape in, 01875 KernelIterator kernels, KernelIterator kernelsEnd, 01876 OutIterator outs) 01877 { 01878 vigra_precondition(kernels != kernelsEnd, 01879 "FFTWConvolvePlan::checkShapesComplex(): empty kernel sequence."); 01880 01881 Shape kernelShape = kernels->shape(); 01882 for(; kernels != kernelsEnd; ++kernels, ++outs) 01883 { 01884 vigra_precondition(in == outs->shape(), 01885 "FFTWConvolvePlan::checkShapesComplex(): shape mismatch between input and (one) output."); 01886 if(useFourierKernel) 01887 { 01888 vigra_precondition(kernelShape == kernels->shape(), 01889 "FFTWConvolvePlan::checkShapesComplex(): Fourier domain kernels must have identical size."); 01890 } 01891 else 01892 { 01893 kernelShape = max(kernelShape, kernels->shape()); 01894 } 01895 } 01896 vigra_precondition(prod(kernelShape) > 0, 01897 "FFTWConvolvePlan::checkShapesComplex(): all kernels have size 0."); 01898 01899 if(useFourierKernel) 01900 { 01901 for(unsigned int k=0; k<N; ++k) 01902 vigra_precondition(in[k] <= kernelShape[k], 01903 "FFTWConvolvePlan::checkShapesComplex(): kernels too small for given input."); 01904 return kernelShape; 01905 } 01906 else 01907 { 01908 return fftwBestPaddedShape(in + kernelShape - Shape(1)); 01909 } 01910 } 01911 01912 01913 /********************************************************/ 01914 /* */ 01915 /* fourierTransform */ 01916 /* */ 01917 /********************************************************/ 01918 01919 template <unsigned int N, class Real, class C1, class C2> 01920 inline void 01921 fourierTransform(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01922 MultiArrayView<N, FFTWComplex<Real>, C2> out) 01923 { 01924 FFTWPlan<N, Real>(in, out, FFTW_FORWARD).execute(in, out); 01925 } 01926 01927 template <unsigned int N, class Real, class C1, class C2> 01928 inline void 01929 fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01930 MultiArrayView<N, FFTWComplex<Real>, C2> out) 01931 { 01932 FFTWPlan<N, Real>(in, out, FFTW_BACKWARD).execute(in, out); 01933 } 01934 01935 template <unsigned int N, class Real, class C1, class C2> 01936 void 01937 fourierTransform(MultiArrayView<N, Real, C1> in, 01938 MultiArrayView<N, FFTWComplex<Real>, C2> out) 01939 { 01940 if(in.shape() == out.shape()) 01941 { 01942 // copy the input array into the output and then perform an in-place FFT 01943 out = in; 01944 FFTWPlan<N, Real>(out, out, FFTW_FORWARD).execute(out, out); 01945 } 01946 else if(out.shape() == fftwCorrespondingShapeR2C(in.shape())) 01947 { 01948 FFTWPlan<N, Real>(in, out).execute(in, out); 01949 } 01950 else 01951 vigra_precondition(false, 01952 "fourierTransform(): shape mismatch between input and output."); 01953 } 01954 01955 template <unsigned int N, class Real, class C1, class C2> 01956 void 01957 fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01958 MultiArrayView<N, Real, C2> out) 01959 { 01960 vigra_precondition(in.shape() == fftwCorrespondingShapeR2C(out.shape()), 01961 "fourierTransformInverse(): shape mismatch between input and output."); 01962 FFTWPlan<N, Real>(in, out).execute(in, out); 01963 } 01964 01965 //@} 01966 01967 /** \addtogroup MultiArrayConvolutionFilters 01968 */ 01969 //@{ 01970 01971 /********************************************************/ 01972 /* */ 01973 /* convolveFFT */ 01974 /* */ 01975 /********************************************************/ 01976 01977 /** \brief Convolve an array with a kernel by means of the Fourier transform. 01978 01979 Thanks to the convolution theorem of Fourier theory, a convolution in the spatial domain 01980 is equivalent to a multiplication in the frequency domain. Thus, for certain kernels 01981 (especially large, non-separable ones), it is advantageous to perform the convolution by first 01982 transforming both array and kernel to the frequency domain, multiplying the frequency 01983 representations, and transforming the result back into the spatial domain. 01984 Some kernels have a much simpler definition in the frequency domain, so that they are readily 01985 computed there directly, avoiding Fourier transformation of those kernels. 01986 01987 The following functions implement various variants of FFT-based convolution: 01988 01989 <DL> 01990 <DT><b>convolveFFT</b><DD> Convolve a real-valued input array with a kernel such that the 01991 result is also real-valued. That is, the kernel is either provided 01992 as a real-valued array in the spatial domain, or as a 01993 complex-valued array in the Fourier domain, using the half-space format 01994 of the R2C Fourier transform (see below). 01995 <DT><b>convolveFFTMany</b><DD> Like <tt>convolveFFT</tt>, but you may provide many kernels at once 01996 (using an iterator pair specifying the kernel sequence). 01997 This has the advantage that the forward transform of the input array needs 01998 to be executed only once. 01999 <DT><b>convolveFFTComplex</b><DD> Convolve a complex-valued input array with a complex-valued kernel, 02000 resulting in a complex-valued output array. An additional flag is used to 02001 specify whether the kernel is defined in the spatial or frequency domain. 02002 <DT><b>convolveFFTComplexMany</b><DD> Like <tt>convolveFFTComplex</tt>, but you may provide many kernels at once 02003 (using an iterator pair specifying the kernel sequence). 02004 This has the advantage that the forward transform of the input array needs 02005 to be executed only once. 02006 </DL> 02007 02008 The output arrays must have the same shape as the input arrays. In the "Many" variants of the 02009 convolution functions, the kernels must all have the same shape. 02010 02011 The origin of the kernel is always assumed to be in the center of the kernel array (precisely, 02012 at the point <tt>floor(kernel.shape() / 2.0)</tt>, except when the half-space format is used, see below). 02013 The function \ref moveDCToUpperLeft() will be called internally to align the kernel with the transformed 02014 input as appropriate. 02015 02016 If a real input is combined with a real kernel, the kernel is automatically assumed to be defined 02017 in the spatial domain. If a real input is combined with a complex kernel, the kernel is assumed 02018 to be defined in the Fourier domain in half-space format. If the input array is complex, a flag 02019 <tt>fourierDomainKernel</tt> determines where the kernel is defined. 02020 02021 When the kernel is defined in the spatial domain, the convolution functions will automatically pad 02022 (enlarge) the input array by at least the kernel radius in each direction. The newly added space is 02023 filled according to reflective boundary conditions in order to minimize border artifacts during 02024 convolution. It is thus ensured that convolution in the Fourier domain yields the same results as 02025 convolution in the spatial domain (e.g. when \ref separableConvolveMultiArray() is called with the 02026 same kernel). A little further padding may be added to make sure that the padded array shape 02027 uses integers which have only small prime factors, because FFTW is then able to use the fastest 02028 possible algorithms. Any padding is automatically removed from the result arrays before the function 02029 returns. 02030 02031 When the kernel is defined in the frequency domain, it must be complex-valued, and its shape determines 02032 the shape of the Fourier representation (i.e. the input is padded according to the shape of the kernel). 02033 If we are going to perform a complex-valued convolution, the kernel must be defined for the entire 02034 frequency domain, and its shape directly determines the size of the FFT. 02035 02036 In contrast, a frequency domain kernel for a real-valued convolution must have symmetry properties 02037 that allow to drop half of the kernel coefficients, as in the 02038 <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C transform</a>. 02039 That is, the kernel must have the <i>half-space format</i>, that is the shape returned by <tt>fftwCorrespondingShapeR2C(fourier_shape)</tt>, where <tt>fourier_shape</tt> is the desired 02040 logical shape of the frequency representation (and thus the size of the padded input). The origin 02041 of the kernel must be at the point 02042 <tt>(0, floor(fourier_shape[0] / 2.0), ..., floor(fourier_shape[N-1] / 2.0))</tt> 02043 (i.e. as in a regular kernel except for the first dimension). 02044 02045 The <tt>Real</tt> type in the declarations can be <tt>double</tt>, <tt>float</tt>, and 02046 <tt>long double</tt>. Your program must always link against <tt>libfftw3</tt>. If you use 02047 <tt>float</tt> or <tt>long double</tt> arrays, you must <i>additionally</i> link against 02048 <tt>libfftw3f</tt> and <tt>libfftw3l</tt> respectively. 02049 02050 The Fourier transform functions internally create <a href="http://www.fftw.org/doc/Using-Plans.html">FFTW plans</a> 02051 which control the algorithm details. The plans are creates with the flag <tt>FFTW_ESTIMATE</tt>, i.e. 02052 optimal settings are guessed or read from saved "wisdom" files. If you need more control over planning, 02053 you can use the class \ref FFTWConvolvePlan. 02054 02055 See also \ref applyFourierFilter() for corresponding functionality on the basis of the 02056 old image iterator interface. 02057 02058 <b> Declarations:</b> 02059 02060 Real-valued convolution with kernel in the spatial domain: 02061 \code 02062 namespace vigra { 02063 template <unsigned int N, class Real, class C1, class C2, class C3> 02064 void 02065 convolveFFT(MultiArrayView<N, Real, C1> in, 02066 MultiArrayView<N, Real, C2> kernel, 02067 MultiArrayView<N, Real, C3> out); 02068 } 02069 \endcode 02070 02071 Real-valued convolution with kernel in the Fourier domain (half-space format): 02072 \code 02073 namespace vigra { 02074 template <unsigned int N, class Real, class C1, class C2, class C3> 02075 void 02076 convolveFFT(MultiArrayView<N, Real, C1> in, 02077 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 02078 MultiArrayView<N, Real, C3> out); 02079 } 02080 \endcode 02081 02082 Series of real-valued convolutions with kernels in the spatial or Fourier domain 02083 (the kernel and out sequences must have the same length): 02084 \code 02085 namespace vigra { 02086 template <unsigned int N, class Real, class C1, 02087 class KernelIterator, class OutIterator> 02088 void 02089 convolveFFTMany(MultiArrayView<N, Real, C1> in, 02090 KernelIterator kernels, KernelIterator kernelsEnd, 02091 OutIterator outs); 02092 } 02093 \endcode 02094 02095 Complex-valued convolution (parameter <tt>fourierDomainKernel</tt> determines if 02096 the kernel is defined in the spatial or Fourier domain): 02097 \code 02098 namespace vigra { 02099 template <unsigned int N, class Real, class C1, class C2, class C3> 02100 void 02101 convolveFFTComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in, 02102 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 02103 MultiArrayView<N, FFTWComplex<Real>, C3> out, 02104 bool fourierDomainKernel); 02105 } 02106 \endcode 02107 02108 Series of complex-valued convolutions (parameter <tt>fourierDomainKernel</tt> 02109 determines if the kernels are defined in the spatial or Fourier domain, 02110 the kernel and out sequences must have the same length): 02111 \code 02112 namespace vigra { 02113 template <unsigned int N, class Real, class C1, 02114 class KernelIterator, class OutIterator> 02115 void 02116 convolveFFTComplexMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 02117 KernelIterator kernels, KernelIterator kernelsEnd, 02118 OutIterator outs, 02119 bool fourierDomainKernel); 02120 } 02121 \endcode 02122 02123 <b> Usage:</b> 02124 02125 <b>\#include</b> <vigra/multi_fft.hxx><br> 02126 Namespace: vigra 02127 02128 \code 02129 // convolve real array with a Gaussian (sigma=1) defined in the spatial domain 02130 // (implicitly uses padding by at least 4 pixels) 02131 MultiArray<2, double> src(Shape2(w, h)), dest(Shape2(w,h)); 02132 02133 MultiArray<2, double> spatial_kernel(Shape2(9, 9)); 02134 Gaussian<double> gauss(1.0); 02135 02136 for(int y=0; y<9; ++y) 02137 for(int x=0; x<9; ++x) 02138 spatial_kernel(x, y) = gauss(x-4.0)*gauss(y-4.0); 02139 02140 convolveFFT(src, spatial_kernel, dest); 02141 02142 // convolve real array with a Gaussian (sigma=1) defined in the Fourier domain 02143 // (uses no padding, because the kernel size corresponds to the input size) 02144 MultiArray<2, FFTWComplex<double> > fourier_kernel(fftwCorrespondingShapeR2C(src.shape())); 02145 int y0 = h / 2; 02146 02147 for(int y=0; y<fourier_kernel.shape(1); ++y) 02148 for(int x=0; x<fourier_kernel.shape(0); ++x) 02149 fourier_kernel(x, y) = exp(-0.5*sq(x / double(w))) * exp(-0.5*sq((y-y0)/double(h))); 02150 02151 convolveFFT(src, fourier_kernel, dest); 02152 \endcode 02153 */ 02154 doxygen_overloaded_function(template <...> void convolveFFT) 02155 02156 template <unsigned int N, class Real, class C1, class C2, class C3> 02157 void 02158 convolveFFT(MultiArrayView<N, Real, C1> in, 02159 MultiArrayView<N, Real, C2> kernel, 02160 MultiArrayView<N, Real, C3> out) 02161 { 02162 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out); 02163 } 02164 02165 template <unsigned int N, class Real, class C1, class C2, class C3> 02166 void 02167 convolveFFT(MultiArrayView<N, Real, C1> in, 02168 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 02169 MultiArrayView<N, Real, C3> out) 02170 { 02171 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out); 02172 } 02173 02174 /** \brief Convolve a complex-valued array by means of the Fourier transform. 02175 02176 See \ref convolveFFT() for details. 02177 */ 02178 doxygen_overloaded_function(template <...> void convolveFFTComplex) 02179 02180 template <unsigned int N, class Real, class C1, class C2, class C3> 02181 void 02182 convolveFFTComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in, 02183 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 02184 MultiArrayView<N, FFTWComplex<Real>, C3> out, 02185 bool fourierDomainKernel) 02186 { 02187 FFTWConvolvePlan<N, Real>(in, kernel, out, fourierDomainKernel).execute(in, kernel, out); 02188 } 02189 02190 /** \brief Convolve a real-valued array with a sequence of kernels by means of the Fourier transform. 02191 02192 See \ref convolveFFT() for details. 02193 */ 02194 doxygen_overloaded_function(template <...> void convolveFFTMany) 02195 02196 template <unsigned int N, class Real, class C1, 02197 class KernelIterator, class OutIterator> 02198 void 02199 convolveFFTMany(MultiArrayView<N, Real, C1> in, 02200 KernelIterator kernels, KernelIterator kernelsEnd, 02201 OutIterator outs) 02202 { 02203 FFTWConvolvePlan<N, Real> plan; 02204 plan.initMany(in, kernels, kernelsEnd, outs); 02205 plan.executeMany(in, kernels, kernelsEnd, outs); 02206 } 02207 02208 /** \brief Convolve a complex-valued array with a sequence of kernels by means of the Fourier transform. 02209 02210 See \ref convolveFFT() for details. 02211 */ 02212 doxygen_overloaded_function(template <...> void convolveFFTComplexMany) 02213 02214 template <unsigned int N, class Real, class C1, 02215 class KernelIterator, class OutIterator> 02216 void 02217 convolveFFTComplexMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 02218 KernelIterator kernels, KernelIterator kernelsEnd, 02219 OutIterator outs, 02220 bool fourierDomainKernel) 02221 { 02222 FFTWConvolvePlan<N, Real> plan; 02223 plan.initMany(in, kernels, kernelsEnd, outs, fourierDomainKernel); 02224 plan.executeMany(in, kernels, kernelsEnd, outs); 02225 } 02226 02227 //@} 02228 02229 } // namespace vigra 02230 02231 #endif // VIGRA_MULTI_FFT_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|