[ 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 // find the smallest padding size that is at least as large as 's' 00499 int * upperBound = std::upper_bound(goodSizes, goodSizes+size, s); 00500 if(upperBound > goodSizes && upperBound[-1] == s) 00501 return s; 00502 else 00503 return *upperBound; 00504 } 00505 00506 inline 00507 int fftwEvenPaddingSize(int s) 00508 { 00509 // Image sizes where FFTW is fast. The list contains all even 00510 // numbers less than 100000 whose prime decomposition is of the form 00511 // 2^a*3^b*5^c*7^d*11^e*13^f, where a >= 1, e+f is either 0 or 1, and the 00512 // other exponents are arbitrary 00513 static const int size = 1063; 00514 static int goodSizes[size] = { 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 00515 24, 26, 28, 30, 32, 36, 40, 42, 44, 48, 50, 52, 54, 56, 60, 64, 66, 70, 00516 72, 78, 80, 84, 88, 90, 96, 98, 100, 104, 108, 110, 112, 120, 126, 128, 00517 130, 132, 140, 144, 150, 154, 156, 160, 162, 168, 176, 180, 182, 192, 00518 196, 198, 200, 208, 210, 216, 220, 224, 234, 240, 250, 252, 256, 260, 00519 264, 270, 280, 288, 294, 300, 308, 312, 320, 324, 330, 336, 350, 352, 00520 360, 364, 378, 384, 390, 392, 396, 400, 416, 420, 432, 440, 448, 450, 00521 462, 468, 480, 486, 490, 500, 504, 512, 520, 528, 540, 546, 550, 560, 00522 576, 588, 594, 600, 616, 624, 630, 640, 648, 650, 660, 672, 686, 700, 00523 702, 704, 720, 728, 750, 756, 768, 770, 780, 784, 792, 800, 810, 832, 00524 840, 864, 880, 882, 896, 900, 910, 924, 936, 960, 972, 980, 990, 1000, 00525 1008, 1024, 1040, 1050, 1056, 1078, 1080, 1092, 1100, 1120, 1134, 1152, 00526 1170, 1176, 1188, 1200, 1232, 1248, 1250, 1260, 1274, 1280, 1296, 1300, 00527 1320, 1344, 1350, 1372, 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470, 00528 1500, 1512, 1536, 1540, 1560, 1568, 1584, 1600, 1620, 1638, 1650, 1664, 00529 1680, 1728, 1750, 1760, 1764, 1782, 1792, 1800, 1820, 1848, 1872, 1890, 00530 1920, 1944, 1950, 1960, 1980, 2000, 2016, 2048, 2058, 2080, 2100, 2106, 00531 2112, 2156, 2160, 2184, 2200, 2240, 2250, 2268, 2304, 2310, 2340, 2352, 00532 2376, 2400, 2430, 2450, 2464, 2496, 2500, 2520, 2548, 2560, 2592, 2600, 00533 2640, 2646, 2688, 2700, 2730, 2744, 2750, 2772, 2800, 2808, 2816, 2880, 00534 2912, 2916, 2940, 2970, 3000, 3024, 3072, 3080, 3120, 3136, 3150, 3168, 00535 3200, 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3402, 3430, 3456, 3500, 00536 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3696, 3744, 3750, 3780, 3822, 00537 3840, 3850, 3888, 3900, 3920, 3960, 4000, 4032, 4050, 4096, 4116, 4158, 00538 4160, 4200, 4212, 4224, 4312, 4320, 4368, 4374, 4400, 4410, 4480, 4500, 00539 4536, 4550, 4608, 4620, 4680, 4704, 4752, 4800, 4802, 4860, 4900, 4914, 00540 4928, 4950, 4992, 5000, 5040, 5096, 5120, 5184, 5200, 5250, 5280, 5292, 00541 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544, 5600, 5616, 5632, 5670, 00542 5760, 5824, 5832, 5850, 5880, 5940, 6000, 6048, 6144, 6160, 6174, 6240, 00543 6250, 6272, 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6600, 00544 6656, 6720, 6750, 6804, 6860, 6912, 6930, 7000, 7020, 7040, 7056, 7128, 00545 7168, 7200, 7280, 7290, 7350, 7392, 7488, 7500, 7546, 7560, 7644, 7680, 00546 7700, 7776, 7800, 7840, 7920, 7938, 8000, 8064, 8100, 8190, 8192, 8232, 00547 8250, 8316, 8320, 8400, 8424, 8448, 8624, 8640, 8736, 8748, 8750, 8800, 00548 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9360, 9408, 9450, 00549 9504, 9600, 9604, 9702, 9720, 9750, 9800, 9828, 9856, 9900, 9984, 10000, 00550 10080, 10192, 10206, 10240, 10290, 10368, 10400, 10500, 10530, 10560, 00551 10584, 10692, 10752, 10780, 10800, 10920, 10976, 11000, 11088, 11200, 00552 11232, 11250, 11264, 11340, 11466, 11520, 11550, 11648, 11664, 11700, 00553 11760, 11880, 12000, 12096, 12150, 12250, 12288, 12320, 12348, 12474, 00554 12480, 12500, 12544, 12600, 12636, 12672, 12740, 12800, 12936, 12960, 00555 13000, 13104, 13122, 13200, 13230, 13312, 13440, 13500, 13608, 13650, 00556 13720, 13750, 13824, 13860, 14000, 14040, 14080, 14112, 14256, 14336, 00557 14400, 14406, 14560, 14580, 14700, 14742, 14784, 14850, 14976, 15000, 00558 15092, 15120, 15288, 15360, 15400, 15552, 15600, 15680, 15750, 15840, 00559 15876, 16000, 16038, 16128, 16170, 16200, 16250, 16380, 16384, 16464, 00560 16500, 16632, 16640, 16800, 16848, 16896, 17010, 17150, 17248, 17280, 00561 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000, 00562 18144, 18200, 18432, 18480, 18522, 18720, 18750, 18816, 18900, 18954, 00563 19008, 19110, 19200, 19208, 19250, 19404, 19440, 19500, 19600, 19656, 00564 19712, 19800, 19968, 20000, 20160, 20250, 20384, 20412, 20480, 20580, 00565 20736, 20790, 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560, 00566 21600, 21840, 21870, 21952, 22000, 22050, 22176, 22400, 22464, 22500, 00567 22528, 22638, 22680, 22750, 22932, 23040, 23100, 23296, 23328, 23400, 00568 23520, 23760, 23814, 24000, 24010, 24192, 24300, 24500, 24570, 24576, 00569 24640, 24696, 24750, 24948, 24960, 25000, 25088, 25200, 25272, 25344, 00570 25480, 25600, 25872, 25920, 26000, 26208, 26244, 26250, 26400, 26460, 00571 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440, 27500, 00572 27648, 27720, 28000, 28080, 28160, 28224, 28350, 28512, 28672, 28800, 00573 28812, 29106, 29120, 29160, 29250, 29400, 29484, 29568, 29700, 29952, 00574 30000, 30184, 30240, 30576, 30618, 30720, 30800, 30870, 31104, 31200, 00575 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256, 00576 32340, 32400, 32500, 32760, 32768, 32928, 33000, 33264, 33280, 33600, 00577 33614, 33696, 33750, 33792, 34020, 34300, 34398, 34496, 34560, 34650, 00578 34944, 34992, 35000, 35100, 35200, 35280, 35640, 35672, 35840, 36000, 00579 36288, 36400, 36450, 36750, 36864, 36960, 37044, 37422, 37440, 37500, 00580 37632, 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808, 00581 38880, 39000, 39200, 39312, 39366, 39424, 39600, 39690, 39936, 40000, 00582 40320, 40500, 40768, 40824, 40950, 40960, 41160, 41250, 41472, 41580, 00583 41600, 42000, 42120, 42240, 42336, 42768, 43008, 43120, 43200, 43218, 00584 43680, 43740, 43750, 43904, 44000, 44100, 44226, 44352, 44550, 44590, 00585 44800, 44928, 45000, 45056, 45276, 45360, 45500, 45864, 46080, 46200, 00586 46592, 46656, 46800, 47040, 47250, 47520, 47628, 48000, 48020, 48114, 00587 48384, 48510, 48600, 48750, 49000, 49140, 49152, 49280, 49392, 49500, 00588 49896, 49920, 50000, 50176, 50400, 50544, 50688, 50960, 51030, 51200, 00589 51450, 51744, 51840, 52000, 52416, 52488, 52500, 52650, 52800, 52822, 00590 52920, 53248, 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54880, 00591 55000, 55296, 55440, 55566, 56000, 56160, 56250, 56320, 56448, 56700, 00592 56862, 57024, 57330, 57344, 57600, 57624, 57750, 58212, 58240, 58320, 00593 58500, 58800, 58968, 59136, 59400, 59904, 60000, 60368, 60480, 60750, 00594 61152, 61236, 61250, 61440, 61600, 61740, 62208, 62370, 62400, 62426, 00595 62500, 62720, 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512, 00596 64680, 64800, 65000, 65520, 65536, 65610, 65856, 66000, 66150, 66528, 00597 66560, 67200, 67228, 67392, 67500, 67584, 67914, 68040, 68250, 68600, 00598 68750, 68796, 68992, 69120, 69300, 69888, 69984, 70000, 70200, 70400, 00599 70560, 71280, 71344, 71442, 71680, 72000, 72030, 72576, 72800, 72900, 00600 73500, 73710, 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264, 00601 75460, 75600, 75816, 76032, 76440, 76800, 76832, 77000, 77616, 77760, 00602 78000, 78400, 78624, 78732, 78750, 78848, 79200, 79380, 79872, 80000, 00603 80190, 80262, 80640, 80850, 81000, 81250, 81536, 81648, 81900, 81920, 00604 82320, 82500, 82944, 83160, 83200, 84000, 84240, 84480, 84672, 85050, 00605 85536, 85750, 86016, 86240, 86400, 86436, 87318, 87360, 87480, 87500, 00606 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856, 00607 90000, 90112, 90552, 90720, 91000, 91728, 91854, 92160, 92400, 92610, 00608 93184, 93312, 93600, 93750, 94080, 94500, 94770, 95040, 95256, 95550, 00609 96000, 96040, 96228, 96250, 96768, 97020, 97200, 97500, 98000, 98280, 00610 98304, 98560, 98784, 99000, 99792, 99840 }; 00611 00612 if(s <= 0 || s >= goodSizes[size-1]) 00613 return s; 00614 // find the smallest padding size that is at least as large as 's' 00615 int * upperBound = std::upper_bound(goodSizes, goodSizes+size, s); 00616 if(upperBound > goodSizes && upperBound[-1] == s) 00617 return s; 00618 else 00619 return *upperBound; 00620 } 00621 00622 template <int M> 00623 struct FFTEmbedKernel 00624 { 00625 template <unsigned int N, class Real, class C, class Shape> 00626 static void 00627 exec(MultiArrayView<N, Real, C> & out, Shape const & kernelShape, 00628 Shape & srcPoint, Shape & destPoint, bool copyIt) 00629 { 00630 for(srcPoint[M]=0; srcPoint[M]<kernelShape[M]; ++srcPoint[M]) 00631 { 00632 if(srcPoint[M] < (kernelShape[M] + 1) / 2) 00633 { 00634 destPoint[M] = srcPoint[M]; 00635 } 00636 else 00637 { 00638 destPoint[M] = srcPoint[M] + out.shape(M) - kernelShape[M]; 00639 copyIt = true; 00640 } 00641 FFTEmbedKernel<M-1>::exec(out, kernelShape, srcPoint, destPoint, copyIt); 00642 } 00643 } 00644 }; 00645 00646 template <> 00647 struct FFTEmbedKernel<0> 00648 { 00649 template <unsigned int N, class Real, class C, class Shape> 00650 static void 00651 exec(MultiArrayView<N, Real, C> & out, Shape const & kernelShape, 00652 Shape & srcPoint, Shape & destPoint, bool copyIt) 00653 { 00654 for(srcPoint[0]=0; srcPoint[0]<kernelShape[0]; ++srcPoint[0]) 00655 { 00656 if(srcPoint[0] < (kernelShape[0] + 1) / 2) 00657 { 00658 destPoint[0] = srcPoint[0]; 00659 } 00660 else 00661 { 00662 destPoint[0] = srcPoint[0] + out.shape(0) - kernelShape[0]; 00663 copyIt = true; 00664 } 00665 if(copyIt) 00666 { 00667 out[destPoint] = out[srcPoint]; 00668 out[srcPoint] = 0.0; 00669 } 00670 } 00671 } 00672 }; 00673 00674 template <unsigned int N, class Real, class C1, class C2> 00675 void 00676 fftEmbedKernel(MultiArrayView<N, Real, C1> kernel, 00677 MultiArrayView<N, Real, C2> out, 00678 Real norm = 1.0) 00679 { 00680 typedef typename MultiArrayShape<N>::type Shape; 00681 00682 MultiArrayView<N, Real, C2> kout = out.subarray(Shape(), kernel.shape()); 00683 00684 out.init(0.0); 00685 kout = kernel; 00686 if (norm != 1.0) 00687 kout *= norm; 00688 moveDCToUpperLeft(kout); 00689 00690 Shape srcPoint, destPoint; 00691 FFTEmbedKernel<(int)N-1>::exec(out, kernel.shape(), srcPoint, destPoint, false); 00692 } 00693 00694 template <unsigned int N, class Real, class C1, class C2> 00695 void 00696 fftEmbedArray(MultiArrayView<N, Real, C1> in, 00697 MultiArrayView<N, Real, C2> out) 00698 { 00699 typedef typename MultiArrayShape<N>::type Shape; 00700 00701 Shape diff = out.shape() - in.shape(), 00702 leftDiff = div(diff, MultiArrayIndex(2)), 00703 rightDiff = diff - leftDiff, 00704 right = in.shape() + leftDiff; 00705 00706 out.subarray(leftDiff, right) = in; 00707 00708 typedef typename MultiArrayView<N, Real, C2>::traverser Traverser; 00709 typedef MultiArrayNavigator<Traverser, N> Navigator; 00710 typedef typename Navigator::iterator Iterator; 00711 00712 for(unsigned int d = 0; d < N; ++d) 00713 { 00714 Navigator nav(out.traverser_begin(), out.shape(), d); 00715 00716 for( ; nav.hasMore(); nav++ ) 00717 { 00718 Iterator i = nav.begin(); 00719 for(int k=1; k<=leftDiff[d]; ++k) 00720 i[leftDiff[d] - k] = i[leftDiff[d] + k]; 00721 for(int k=0; k<rightDiff[d]; ++k) 00722 i[right[d] + k] = i[right[d] - k - 2]; 00723 } 00724 } 00725 } 00726 00727 } // namespace detail 00728 00729 template <class T, int N> 00730 TinyVector<T, N> 00731 fftwBestPaddedShape(TinyVector<T, N> shape) 00732 { 00733 for(unsigned int k=0; k<N; ++k) 00734 shape[k] = detail::fftwPaddingSize(shape[k]); 00735 return shape; 00736 } 00737 00738 template <class T, int N> 00739 TinyVector<T, N> 00740 fftwBestPaddedShapeR2C(TinyVector<T, N> shape) 00741 { 00742 shape[0] = detail::fftwEvenPaddingSize(shape[0]); 00743 for(unsigned int k=1; k<N; ++k) 00744 shape[k] = detail::fftwPaddingSize(shape[k]); 00745 return shape; 00746 } 00747 00748 /** \brief Find frequency domain shape for a R2C Fourier transform. 00749 00750 When a real valued array is transformed to the frequency domain, about half of the 00751 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 00752 transform</a> that doesn't compute and store the redundant coefficients. This function 00753 computes the appropriate frequency domain shape for a given shape in the spatial domain. 00754 It simply replaces <tt>shape[0]</tt> with <tt>shape[0] / 2 + 1</tt>. 00755 00756 <b>\#include</b> <vigra/multi_fft.hxx> 00757 */ 00758 template <class T, int N> 00759 TinyVector<T, N> 00760 fftwCorrespondingShapeR2C(TinyVector<T, N> shape) 00761 { 00762 shape[0] = shape[0] / 2 + 1; 00763 return shape; 00764 } 00765 00766 template <class T, int N> 00767 TinyVector<T, N> 00768 fftwCorrespondingShapeC2R(TinyVector<T, N> shape, bool oddDimension0 = false) 00769 { 00770 shape[0] = oddDimension0 00771 ? (shape[0] - 1) * 2 + 1 00772 : (shape[0] - 1) * 2; 00773 return shape; 00774 } 00775 00776 /********************************************************/ 00777 /* */ 00778 /* FFTWPlan */ 00779 /* */ 00780 /********************************************************/ 00781 00782 /** C++ wrapper for FFTW plans. 00783 00784 The class encapsulates the calls to <tt>fftw_plan_dft_2d</tt>, <tt>fftw_execute</tt>, and 00785 <tt>fftw_destroy_plan</tt> (and their <tt>float</tt> and <tt>long double</tt> counterparts) 00786 in an easy-to-use interface. 00787 00788 Usually, you use this class only indirectly via \ref fourierTransform() 00789 and \ref fourierTransformInverse(). You only need this class if you want to have more control 00790 about FFTW's planning process (by providing non-default planning flags) and/or want to re-use 00791 plans for several transformations. 00792 00793 <b> Usage:</b> 00794 00795 <b>\#include</b> <vigra/multi_fft.hxx><br> 00796 Namespace: vigra 00797 00798 \code 00799 // compute complex Fourier transform of a real image 00800 MultiArray<2, double> src(Shape2(w, h)); 00801 MultiArray<2, FFTWComplex<double> > fourier(Shape2(w, h)); 00802 00803 // create an optimized plan by measuring the speed of several algorithm variants 00804 FFTWPlan<2, double> plan(src, fourier, FFTW_MEASURE); 00805 00806 plan.execute(src, fourier); 00807 \endcode 00808 */ 00809 template <unsigned int N, class Real = double> 00810 class FFTWPlan 00811 { 00812 typedef ArrayVector<int> Shape; 00813 typedef typename FFTWReal2Complex<Real>::plan_type PlanType; 00814 typedef typename FFTWComplex<Real>::complex_type Complex; 00815 00816 PlanType plan; 00817 Shape shape, instrides, outstrides; 00818 int sign; 00819 00820 public: 00821 /** \brief Create an empty plan. 00822 00823 The plan can be initialized later by one of the init() functions. 00824 */ 00825 FFTWPlan() 00826 : plan(0) 00827 {} 00828 00829 /** \brief Create a plan for a complex-to-complex transform. 00830 00831 \arg SIGN must be <tt>FFTW_FORWARD</tt> or <tt>FFTW_BACKWARD</tt> according to the 00832 desired transformation direction. 00833 \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 00834 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 00835 optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 00836 */ 00837 template <class C1, class C2> 00838 FFTWPlan(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00839 MultiArrayView<N, FFTWComplex<Real>, C2> out, 00840 int SIGN, unsigned int planner_flags = FFTW_ESTIMATE) 00841 : plan(0) 00842 { 00843 init(in, out, SIGN, planner_flags); 00844 } 00845 00846 /** \brief Create a plan for a real-to-complex transform. 00847 00848 This always refers to a forward transform. The shape of the output determines 00849 if a standard transform (when <tt>out.shape() == in.shape()</tt>) or an 00850 <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C 00851 transform</a> (when <tt>out.shape() == fftwCorrespondingShapeR2C(in.shape())</tt>) will be executed. 00852 00853 \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 00854 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 00855 optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 00856 */ 00857 template <class C1, class C2> 00858 FFTWPlan(MultiArrayView<N, Real, C1> in, 00859 MultiArrayView<N, FFTWComplex<Real>, C2> out, 00860 unsigned int planner_flags = FFTW_ESTIMATE) 00861 : plan(0) 00862 { 00863 init(in, out, planner_flags); 00864 } 00865 00866 /** \brief Create a plan for a complex-to-real transform. 00867 00868 This always refers to a inverse transform. The shape of the input determines 00869 if a standard transform (when <tt>in.shape() == out.shape()</tt>) or a 00870 <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">C2R 00871 transform</a> (when <tt>in.shape() == fftwCorrespondingShapeR2C(out.shape())</tt>) will be executed. 00872 00873 \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 00874 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 00875 optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 00876 */ 00877 template <class C1, class C2> 00878 FFTWPlan(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00879 MultiArrayView<N, Real, C2> out, 00880 unsigned int planner_flags = FFTW_ESTIMATE) 00881 : plan(0) 00882 { 00883 init(in, out, planner_flags); 00884 } 00885 00886 /** \brief Copy constructor. 00887 */ 00888 FFTWPlan(FFTWPlan const & other) 00889 : plan(other.plan), 00890 sign(other.sign) 00891 { 00892 FFTWPlan & o = const_cast<FFTWPlan &>(other); 00893 shape.swap(o.shape); 00894 instrides.swap(o.instrides); 00895 outstrides.swap(o.outstrides); 00896 o.plan = 0; // act like std::auto_ptr 00897 } 00898 00899 /** \brief Copy assigment. 00900 */ 00901 FFTWPlan & operator=(FFTWPlan const & other) 00902 { 00903 if(this != &other) 00904 { 00905 FFTWPlan & o = const_cast<FFTWPlan &>(other); 00906 plan = o.plan; 00907 shape.swap(o.shape); 00908 instrides.swap(o.instrides); 00909 outstrides.swap(o.outstrides); 00910 sign = o.sign; 00911 o.plan = 0; // act like std::auto_ptr 00912 } 00913 return *this; 00914 } 00915 00916 /** \brief Destructor. 00917 */ 00918 ~FFTWPlan() 00919 { 00920 detail::fftwPlanDestroy(plan); 00921 } 00922 00923 /** \brief Init a complex-to-complex transform. 00924 00925 See the constructor with the same signature for details. 00926 */ 00927 template <class C1, class C2> 00928 void init(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00929 MultiArrayView<N, FFTWComplex<Real>, C2> out, 00930 int SIGN, unsigned int planner_flags = FFTW_ESTIMATE) 00931 { 00932 vigra_precondition(in.strideOrdering() == out.strideOrdering(), 00933 "FFTWPlan.init(): input and output must have the same stride ordering."); 00934 00935 initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(), 00936 SIGN, planner_flags); 00937 } 00938 00939 /** \brief Init a real-to-complex transform. 00940 00941 See the constructor with the same signature for details. 00942 */ 00943 template <class C1, class C2> 00944 void init(MultiArrayView<N, Real, C1> in, 00945 MultiArrayView<N, FFTWComplex<Real>, C2> out, 00946 unsigned int planner_flags = FFTW_ESTIMATE) 00947 { 00948 vigra_precondition(in.strideOrdering() == out.strideOrdering(), 00949 "FFTWPlan.init(): input and output must have the same stride ordering."); 00950 00951 initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(), 00952 FFTW_FORWARD, planner_flags); 00953 } 00954 00955 /** \brief Init a complex-to-real transform. 00956 00957 See the constructor with the same signature for details. 00958 */ 00959 template <class C1, class C2> 00960 void init(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00961 MultiArrayView<N, Real, C2> out, 00962 unsigned int planner_flags = FFTW_ESTIMATE) 00963 { 00964 vigra_precondition(in.strideOrdering() == out.strideOrdering(), 00965 "FFTWPlan.init(): input and output must have the same stride ordering."); 00966 00967 initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(), 00968 FFTW_BACKWARD, planner_flags); 00969 } 00970 00971 /** \brief Execute a complex-to-complex transform. 00972 00973 The array shapes must be the same as in the corresponding init function 00974 or constructor. However, execute() can be called several times on 00975 the same plan, even with different arrays, as long as they have the appropriate 00976 shapes. 00977 */ 00978 template <class C1, class C2> 00979 void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 00980 MultiArrayView<N, FFTWComplex<Real>, C2> out) const 00981 { 00982 executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending()); 00983 } 00984 00985 /** \brief Execute a real-to-complex transform. 00986 00987 The array shapes must be the same as in the corresponding init function 00988 or constructor. However, execute() can be called several times on 00989 the same plan, even with different arrays, as long as they have the appropriate 00990 shapes. 00991 */ 00992 template <class C1, class C2> 00993 void execute(MultiArrayView<N, Real, C1> in, 00994 MultiArrayView<N, FFTWComplex<Real>, C2> out) const 00995 { 00996 executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending()); 00997 } 00998 00999 /** \brief Execute a complex-to-real transform. 01000 01001 The array shapes must be the same as in the corresponding init function 01002 or constructor. However, execute() can be called several times on 01003 the same plan, even with different arrays, as long as they have the appropriate 01004 shapes. 01005 */ 01006 template <class C1, class C2> 01007 void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01008 MultiArrayView<N, Real, C2> out) const 01009 { 01010 executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending()); 01011 } 01012 01013 private: 01014 01015 template <class MI, class MO> 01016 void initImpl(MI ins, MO outs, int SIGN, unsigned int planner_flags); 01017 01018 template <class MI, class MO> 01019 void executeImpl(MI ins, MO outs) const; 01020 01021 void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> in, 01022 MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> out) const 01023 { 01024 vigra_precondition(in.shape() == out.shape(), 01025 "FFTWPlan.init(): input and output must have the same shape."); 01026 } 01027 01028 void checkShapes(MultiArrayView<N, Real, StridedArrayTag> ins, 01029 MultiArrayView<N, FFTWComplex<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(ins.shape(N-1) / 2 + 1 == outs.shape(N-1), 01035 "FFTWPlan.init(): input and output must have matching shapes."); 01036 } 01037 01038 void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> ins, 01039 MultiArrayView<N, Real, StridedArrayTag> outs) const 01040 { 01041 for(int k=0; k<(int)N-1; ++k) 01042 vigra_precondition(ins.shape(k) == outs.shape(k), 01043 "FFTWPlan.init(): input and output must have matching shapes."); 01044 vigra_precondition(outs.shape(N-1) / 2 + 1 == ins.shape(N-1), 01045 "FFTWPlan.init(): input and output must have matching shapes."); 01046 } 01047 }; 01048 01049 template <unsigned int N, class Real> 01050 template <class MI, class MO> 01051 void 01052 FFTWPlan<N, Real>::initImpl(MI ins, MO outs, int SIGN, unsigned int planner_flags) 01053 { 01054 checkShapes(ins, outs); 01055 01056 typename MultiArrayShape<N>::type logicalShape(SIGN == FFTW_FORWARD 01057 ? ins.shape() 01058 : outs.shape()); 01059 01060 Shape newShape(logicalShape.begin(), logicalShape.end()), 01061 newIStrides(ins.stride().begin(), ins.stride().end()), 01062 newOStrides(outs.stride().begin(), outs.stride().end()), 01063 itotal(ins.shape().begin(), ins.shape().end()), 01064 ototal(outs.shape().begin(), outs.shape().end()); 01065 01066 for(unsigned int j=1; j<N; ++j) 01067 { 01068 itotal[j] = ins.stride(j-1) / ins.stride(j); 01069 ototal[j] = outs.stride(j-1) / outs.stride(j); 01070 } 01071 01072 PlanType newPlan = detail::fftwPlanCreate(N, newShape.begin(), 01073 ins.data(), itotal.begin(), ins.stride(N-1), 01074 outs.data(), ototal.begin(), outs.stride(N-1), 01075 SIGN, planner_flags); 01076 detail::fftwPlanDestroy(plan); 01077 plan = newPlan; 01078 shape.swap(newShape); 01079 instrides.swap(newIStrides); 01080 outstrides.swap(newOStrides); 01081 sign = SIGN; 01082 } 01083 01084 template <unsigned int N, class Real> 01085 template <class MI, class MO> 01086 void FFTWPlan<N, Real>::executeImpl(MI ins, MO outs) const 01087 { 01088 vigra_precondition(plan != 0, "FFTWPlan::execute(): plan is NULL."); 01089 01090 typename MultiArrayShape<N>::type lshape(sign == FFTW_FORWARD 01091 ? ins.shape() 01092 : outs.shape()); 01093 01094 vigra_precondition((lshape == TinyVectorView<int, N>(shape.data())), 01095 "FFTWPlan::execute(): shape mismatch between plan and data."); 01096 vigra_precondition((ins.stride() == TinyVectorView<int, N>(instrides.data())), 01097 "FFTWPlan::execute(): strides mismatch between plan and input data."); 01098 vigra_precondition((outs.stride() == TinyVectorView<int, N>(outstrides.data())), 01099 "FFTWPlan::execute(): strides mismatch between plan and output data."); 01100 01101 detail::fftwPlanExecute(plan, ins.data(), outs.data()); 01102 01103 typedef typename MO::value_type V; 01104 if(sign == FFTW_BACKWARD) 01105 outs *= V(1.0) / Real(outs.size()); 01106 } 01107 01108 /********************************************************/ 01109 /* */ 01110 /* FFTWConvolvePlan */ 01111 /* */ 01112 /********************************************************/ 01113 01114 /** C++ wrapper for a pair of FFTW plans used to perform FFT-based convolution. 01115 01116 The class encapsulates the calls to <tt>fftw_plan_dft_2d</tt>, <tt>fftw_execute</tt>, and 01117 <tt>fftw_destroy_plan</tt> (and their <tt>float</tt> and <tt>long double</tt> counterparts) 01118 in an easy-to-use interface. It always creates a pair of plans, one for the forward and one 01119 for the inverse transform required for convolution. 01120 01121 Usually, you use this class only indirectly via \ref convolveFFT() and its variants. 01122 You only need this class if you want to have more control about FFTW's planning process 01123 (by providing non-default planning flags) and/or want to re-use plans for several convolutions. 01124 01125 <b> Usage:</b> 01126 01127 <b>\#include</b> <vigra/multi_fft.hxx><br> 01128 Namespace: vigra 01129 01130 \code 01131 // convolve a real array with a real kernel 01132 MultiArray<2, double> src(Shape2(w, h)), dest(Shape2(w, h)); 01133 01134 MultiArray<2, double> spatial_kernel(Shape2(9, 9)); 01135 Gaussian<double> gauss(1.0); 01136 01137 for(int y=0; y<9; ++y) 01138 for(int x=0; x<9; ++x) 01139 spatial_kernel(x, y) = gauss(x-4.0)*gauss(y-4.0); 01140 01141 // create an optimized plan by measuring the speed of several algorithm variants 01142 FFTWConvolvePlan<2, double> plan(src, spatial_kernel, dest, FFTW_MEASURE); 01143 01144 plan.execute(src, spatial_kernel, dest); 01145 \endcode 01146 */ 01147 template <unsigned int N, class Real> 01148 class FFTWConvolvePlan 01149 { 01150 typedef FFTWComplex<Real> Complex; 01151 typedef MultiArrayView<N, Real, UnstridedArrayTag > RArray; 01152 typedef MultiArray<N, Complex, FFTWAllocator<Complex> > CArray; 01153 01154 FFTWPlan<N, Real> forward_plan, backward_plan; 01155 RArray realArray, realKernel; 01156 CArray fourierArray, fourierKernel; 01157 bool useFourierKernel; 01158 01159 public: 01160 01161 typedef typename MultiArrayShape<N>::type Shape; 01162 01163 /** \brief Create an empty plan. 01164 01165 The plan can be initialized later by one of the init() functions. 01166 */ 01167 FFTWConvolvePlan() 01168 : useFourierKernel(false) 01169 {} 01170 01171 /** \brief Create a plan to convolve a real array with a real kernel. 01172 01173 The kernel must be defined in the spatial domain. 01174 See \ref convolveFFT() for detailed information on required shapes and internal padding. 01175 01176 \arg planner_flags must be a combination of the 01177 <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 01178 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 01179 optimal algorithm settings or read them from pre-loaded 01180 <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 01181 */ 01182 template <class C1, class C2, class C3> 01183 FFTWConvolvePlan(MultiArrayView<N, Real, C1> in, 01184 MultiArrayView<N, Real, C2> kernel, 01185 MultiArrayView<N, Real, C3> out, 01186 unsigned int planner_flags = FFTW_ESTIMATE) 01187 : useFourierKernel(false) 01188 { 01189 init(in, kernel, out, planner_flags); 01190 } 01191 01192 /** \brief Create a plan to convolve a real array with a complex kernel. 01193 01194 The kernel must be defined in the Fourier domain, using the half-space format. 01195 See \ref convolveFFT() for detailed information on required shapes and internal padding. 01196 01197 \arg planner_flags must be a combination of the 01198 <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 01199 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 01200 optimal algorithm settings or read them from pre-loaded 01201 <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 01202 */ 01203 template <class C1, class C2, class C3> 01204 FFTWConvolvePlan(MultiArrayView<N, Real, C1> in, 01205 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01206 MultiArrayView<N, Real, C3> out, 01207 unsigned int planner_flags = FFTW_ESTIMATE) 01208 : useFourierKernel(true) 01209 { 01210 init(in, kernel, out, planner_flags); 01211 } 01212 01213 /** \brief Create a plan to convolve a complex array with a complex kernel. 01214 01215 See \ref convolveFFT() for detailed information on required shapes and internal padding. 01216 01217 \arg fourierDomainKernel determines if the kernel is defined in the spatial or 01218 Fourier domain. 01219 \arg planner_flags must be a combination of the 01220 <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 01221 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 01222 optimal algorithm settings or read them from pre-loaded 01223 <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 01224 */ 01225 template <class C1, class C2, class C3> 01226 FFTWConvolvePlan(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01227 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01228 MultiArrayView<N, FFTWComplex<Real>, C3> out, 01229 bool fourierDomainKernel, 01230 unsigned int planner_flags = FFTW_ESTIMATE) 01231 { 01232 init(in, kernel, out, fourierDomainKernel, planner_flags); 01233 } 01234 01235 01236 /** \brief Create a plan from just the shape information. 01237 01238 See \ref convolveFFT() for detailed information on required shapes and internal padding. 01239 01240 \arg fourierDomainKernel determines if the kernel is defined in the spatial or 01241 Fourier domain. 01242 \arg planner_flags must be a combination of the 01243 <a href="http://www.fftw.org/doc/Planner-Flags.html">planner 01244 flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess 01245 optimal algorithm settings or read them from pre-loaded 01246 <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>. 01247 */ 01248 template <class C1, class C2, class C3> 01249 FFTWConvolvePlan(Shape inOut, Shape kernel, 01250 bool useFourierKernel = false, 01251 unsigned int planner_flags = FFTW_ESTIMATE) 01252 { 01253 if(useFourierKernel) 01254 init(inOut, kernel, planner_flags); 01255 else 01256 initFourierKernel(inOut, kernel, planner_flags); 01257 } 01258 01259 /** \brief Init a plan to convolve a real array with a real kernel. 01260 01261 See the constructor with the same signature for details. 01262 */ 01263 template <class C1, class C2, class C3> 01264 void init(MultiArrayView<N, Real, C1> in, 01265 MultiArrayView<N, Real, C2> kernel, 01266 MultiArrayView<N, Real, C3> out, 01267 unsigned int planner_flags = FFTW_ESTIMATE) 01268 { 01269 vigra_precondition(in.shape() == out.shape(), 01270 "FFTWConvolvePlan::init(): input and output must have the same shape."); 01271 init(in.shape(), kernel.shape(), planner_flags); 01272 } 01273 01274 /** \brief Init a plan to convolve a real array with a complex kernel. 01275 01276 See the constructor with the same signature for details. 01277 */ 01278 template <class C1, class C2, class C3> 01279 void init(MultiArrayView<N, Real, C1> in, 01280 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01281 MultiArrayView<N, Real, C3> out, 01282 unsigned int planner_flags = FFTW_ESTIMATE) 01283 { 01284 vigra_precondition(in.shape() == out.shape(), 01285 "FFTWConvolvePlan::init(): input and output must have the same shape."); 01286 initFourierKernel(in.shape(), kernel.shape(), planner_flags); 01287 } 01288 01289 /** \brief Init a plan to convolve a complex array with a complex kernel. 01290 01291 See the constructor with the same signature for details. 01292 */ 01293 template <class C1, class C2, class C3> 01294 void init(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01295 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01296 MultiArrayView<N, FFTWComplex<Real>, C3> out, 01297 bool fourierDomainKernel, 01298 unsigned int planner_flags = FFTW_ESTIMATE) 01299 { 01300 vigra_precondition(in.shape() == out.shape(), 01301 "FFTWConvolvePlan::init(): input and output must have the same shape."); 01302 useFourierKernel = fourierDomainKernel; 01303 initComplex(in.shape(), kernel.shape(), planner_flags); 01304 } 01305 01306 /** \brief Init a plan to convolve a real array with a sequence of kernels. 01307 01308 The kernels can be either real or complex. The sequences \a kernels and \a outs 01309 must have the same length. See the corresponding constructors 01310 for single kernels for details. 01311 */ 01312 template <class C1, class KernelIterator, class OutIterator> 01313 void initMany(MultiArrayView<N, Real, C1> in, 01314 KernelIterator kernels, KernelIterator kernelsEnd, 01315 OutIterator outs, unsigned int planner_flags = FFTW_ESTIMATE) 01316 { 01317 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray; 01318 typedef typename KernelArray::value_type KernelValue; 01319 typedef typename std::iterator_traits<OutIterator>::value_type OutArray; 01320 typedef typename OutArray::value_type OutValue; 01321 01322 bool realKernel = IsSameType<KernelValue, Real>::value; 01323 bool fourierKernel = IsSameType<KernelValue, Complex>::value; 01324 01325 vigra_precondition(realKernel || fourierKernel, 01326 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type."); 01327 vigra_precondition((IsSameType<OutValue, Real>::value), 01328 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type."); 01329 01330 if(realKernel) 01331 { 01332 initMany(in.shape(), checkShapes(in.shape(), kernels, kernelsEnd, outs), 01333 planner_flags); 01334 } 01335 else 01336 { 01337 initFourierKernelMany(in.shape(), 01338 checkShapesFourier(in.shape(), kernels, kernelsEnd, outs), 01339 planner_flags); 01340 } 01341 } 01342 01343 /** \brief Init a plan to convolve a complex array with a sequence of kernels. 01344 01345 The kernels must be complex as well. The sequences \a kernels and \a outs 01346 must have the same length. See the corresponding constructors 01347 for single kernels for details. 01348 */ 01349 template <class C1, class KernelIterator, class OutIterator> 01350 void initMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01351 KernelIterator kernels, KernelIterator kernelsEnd, 01352 OutIterator outs, 01353 bool fourierDomainKernels, 01354 unsigned int planner_flags = FFTW_ESTIMATE) 01355 { 01356 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray; 01357 typedef typename KernelArray::value_type KernelValue; 01358 typedef typename std::iterator_traits<OutIterator>::value_type OutArray; 01359 typedef typename OutArray::value_type OutValue; 01360 01361 vigra_precondition((IsSameType<KernelValue, Complex>::value), 01362 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type."); 01363 vigra_precondition((IsSameType<OutValue, Complex>::value), 01364 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type."); 01365 01366 useFourierKernel = fourierDomainKernels; 01367 01368 Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs); 01369 01370 CArray newFourierArray(paddedShape), newFourierKernel(paddedShape); 01371 01372 FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags); 01373 FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags); 01374 01375 forward_plan = fplan; 01376 backward_plan = bplan; 01377 fourierArray.swap(newFourierArray); 01378 fourierKernel.swap(newFourierKernel); 01379 } 01380 01381 void init(Shape inOut, Shape kernel, 01382 unsigned int planner_flags = FFTW_ESTIMATE); 01383 01384 void initFourierKernel(Shape inOut, Shape kernel, 01385 unsigned int planner_flags = FFTW_ESTIMATE); 01386 01387 void initComplex(Shape inOut, Shape kernel, 01388 unsigned int planner_flags = FFTW_ESTIMATE); 01389 01390 void initMany(Shape inOut, Shape maxKernel, 01391 unsigned int planner_flags = FFTW_ESTIMATE) 01392 { 01393 init(inOut, maxKernel, planner_flags); 01394 } 01395 01396 void initFourierKernelMany(Shape inOut, Shape kernels, 01397 unsigned int planner_flags = FFTW_ESTIMATE) 01398 { 01399 initFourierKernel(inOut, kernels, planner_flags); 01400 } 01401 01402 /** \brief Execute a plan to convolve a real array with a real kernel. 01403 01404 The array shapes must be the same as in the corresponding init function 01405 or constructor. However, execute() can be called several times on 01406 the same plan, even with different arrays, as long as they have the appropriate 01407 shapes. 01408 */ 01409 template <class C1, class C2, class C3> 01410 void execute(MultiArrayView<N, Real, C1> in, 01411 MultiArrayView<N, Real, C2> kernel, 01412 MultiArrayView<N, Real, C3> out); 01413 01414 /** \brief Execute a plan to convolve a real array with a complex kernel. 01415 01416 The array shapes must be the same as in the corresponding init function 01417 or constructor. However, execute() can be called several times on 01418 the same plan, even with different arrays, as long as they have the appropriate 01419 shapes. 01420 */ 01421 template <class C1, class C2, class C3> 01422 void execute(MultiArrayView<N, Real, C1> in, 01423 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01424 MultiArrayView<N, Real, C3> out); 01425 01426 /** \brief Execute a plan to convolve a complex array with a complex kernel. 01427 01428 The array shapes must be the same as in the corresponding init function 01429 or constructor. However, execute() can be called several times on 01430 the same plan, even with different arrays, as long as they have the appropriate 01431 shapes. 01432 */ 01433 template <class C1, class C2, class C3> 01434 void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01435 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01436 MultiArrayView<N, FFTWComplex<Real>, C3> out); 01437 01438 01439 /** \brief Execute a plan to convolve a complex array with a sequence of kernels. 01440 01441 The array shapes must be the same as in the corresponding init function 01442 or constructor. However, executeMany() can be called several times on 01443 the same plan, even with different arrays, as long as they have the appropriate 01444 shapes. 01445 */ 01446 template <class C1, class KernelIterator, class OutIterator> 01447 void executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01448 KernelIterator kernels, KernelIterator kernelsEnd, 01449 OutIterator outs); 01450 01451 /** \brief Execute a plan to convolve a real array with a sequence of kernels. 01452 01453 The array shapes must be the same as in the corresponding init function 01454 or constructor. However, executeMany() can be called several times on 01455 the same plan, even with different arrays, as long as they have the appropriate 01456 shapes. 01457 */ 01458 template <class C1, class KernelIterator, class OutIterator> 01459 void executeMany(MultiArrayView<N, Real, C1> in, 01460 KernelIterator kernels, KernelIterator kernelsEnd, 01461 OutIterator outs) 01462 { 01463 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray; 01464 typedef typename KernelArray::value_type KernelValue; 01465 typedef typename IsSameType<KernelValue, Complex>::type UseFourierKernel; 01466 typedef typename std::iterator_traits<OutIterator>::value_type OutArray; 01467 typedef typename OutArray::value_type OutValue; 01468 01469 bool realKernel = IsSameType<KernelValue, Real>::value; 01470 bool fourierKernel = IsSameType<KernelValue, Complex>::value; 01471 01472 vigra_precondition(realKernel || fourierKernel, 01473 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type."); 01474 vigra_precondition((IsSameType<OutValue, Real>::value), 01475 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type."); 01476 01477 executeManyImpl(in, kernels, kernelsEnd, outs, UseFourierKernel()); 01478 } 01479 01480 private: 01481 01482 template <class KernelIterator, class OutIterator> 01483 Shape checkShapes(Shape in, 01484 KernelIterator kernels, KernelIterator kernelsEnd, 01485 OutIterator outs); 01486 01487 template <class KernelIterator, class OutIterator> 01488 Shape checkShapesFourier(Shape in, 01489 KernelIterator kernels, KernelIterator kernelsEnd, 01490 OutIterator outs); 01491 01492 template <class KernelIterator, class OutIterator> 01493 Shape checkShapesComplex(Shape in, 01494 KernelIterator kernels, KernelIterator kernelsEnd, 01495 OutIterator outs); 01496 01497 template <class C1, class KernelIterator, class OutIterator> 01498 void 01499 executeManyImpl(MultiArrayView<N, Real, C1> in, 01500 KernelIterator kernels, KernelIterator kernelsEnd, 01501 OutIterator outs, VigraFalseType /* useFourierKernel*/); 01502 01503 template <class C1, class KernelIterator, class OutIterator> 01504 void 01505 executeManyImpl(MultiArrayView<N, Real, C1> in, 01506 KernelIterator kernels, KernelIterator kernelsEnd, 01507 OutIterator outs, VigraTrueType /* useFourierKernel*/); 01508 01509 }; 01510 01511 template <unsigned int N, class Real> 01512 void 01513 FFTWConvolvePlan<N, Real>::init(Shape in, Shape kernel, 01514 unsigned int planner_flags) 01515 { 01516 Shape paddedShape = fftwBestPaddedShapeR2C(in + kernel - Shape(1)), 01517 complexShape = fftwCorrespondingShapeR2C(paddedShape); 01518 01519 CArray newFourierArray(complexShape), newFourierKernel(complexShape); 01520 01521 Shape realStrides = 2*newFourierArray.stride(); 01522 realStrides[0] = 1; 01523 RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data()); 01524 RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data()); 01525 01526 FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags); 01527 FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags); 01528 01529 forward_plan = fplan; 01530 backward_plan = bplan; 01531 realArray = newRealArray; 01532 realKernel = newRealKernel; 01533 fourierArray.swap(newFourierArray); 01534 fourierKernel.swap(newFourierKernel); 01535 useFourierKernel = false; 01536 } 01537 01538 template <unsigned int N, class Real> 01539 void 01540 FFTWConvolvePlan<N, Real>::initFourierKernel(Shape in, Shape kernel, 01541 unsigned int planner_flags) 01542 { 01543 Shape complexShape = kernel, 01544 paddedShape = fftwCorrespondingShapeC2R(complexShape); 01545 01546 for(unsigned int k=0; k<N; ++k) 01547 vigra_precondition(in[k] <= paddedShape[k], 01548 "FFTWConvolvePlan::init(): kernel too small for given input."); 01549 01550 CArray newFourierArray(complexShape), newFourierKernel(complexShape); 01551 01552 Shape realStrides = 2*newFourierArray.stride(); 01553 realStrides[0] = 1; 01554 RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data()); 01555 RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data()); 01556 01557 FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags); 01558 FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags); 01559 01560 forward_plan = fplan; 01561 backward_plan = bplan; 01562 realArray = newRealArray; 01563 realKernel = newRealKernel; 01564 fourierArray.swap(newFourierArray); 01565 fourierKernel.swap(newFourierKernel); 01566 useFourierKernel = true; 01567 } 01568 01569 template <unsigned int N, class Real> 01570 void 01571 FFTWConvolvePlan<N, Real>::initComplex(Shape in, Shape kernel, 01572 unsigned int planner_flags) 01573 { 01574 Shape paddedShape; 01575 01576 if(useFourierKernel) 01577 { 01578 for(unsigned int k=0; k<N; ++k) 01579 vigra_precondition(in[k] <= kernel[k], 01580 "FFTWConvolvePlan::init(): kernel too small for given input."); 01581 01582 paddedShape = kernel; 01583 } 01584 else 01585 { 01586 paddedShape = fftwBestPaddedShape(in + kernel - Shape(1)); 01587 } 01588 01589 CArray newFourierArray(paddedShape), newFourierKernel(paddedShape); 01590 01591 FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags); 01592 FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags); 01593 01594 forward_plan = fplan; 01595 backward_plan = bplan; 01596 fourierArray.swap(newFourierArray); 01597 fourierKernel.swap(newFourierKernel); 01598 } 01599 01600 #ifndef DOXYGEN // doxygen documents these functions as free functions 01601 01602 template <unsigned int N, class Real> 01603 template <class C1, class C2, class C3> 01604 void 01605 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, Real, C1> in, 01606 MultiArrayView<N, Real, C2> kernel, 01607 MultiArrayView<N, Real, C3> out) 01608 { 01609 vigra_precondition(!useFourierKernel, 01610 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel."); 01611 01612 vigra_precondition(in.shape() == out.shape(), 01613 "FFTWConvolvePlan::execute(): input and output must have the same shape."); 01614 01615 Shape paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernel.shape() - Shape(1)), 01616 diff = paddedShape - in.shape(), 01617 left = div(diff, MultiArrayIndex(2)), 01618 right = in.shape() + left; 01619 01620 vigra_precondition(paddedShape == realArray.shape(), 01621 "FFTWConvolvePlan::execute(): shape mismatch between input and plan."); 01622 01623 detail::fftEmbedArray(in, realArray); 01624 forward_plan.execute(realArray, fourierArray); 01625 01626 detail::fftEmbedKernel(kernel, realKernel); 01627 forward_plan.execute(realKernel, fourierKernel); 01628 01629 fourierArray *= fourierKernel; 01630 01631 backward_plan.execute(fourierArray, realArray); 01632 01633 out = realArray.subarray(left, right); 01634 } 01635 01636 template <unsigned int N, class Real> 01637 template <class C1, class C2, class C3> 01638 void 01639 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, Real, C1> in, 01640 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01641 MultiArrayView<N, Real, C3> out) 01642 { 01643 vigra_precondition(useFourierKernel, 01644 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel."); 01645 01646 vigra_precondition(in.shape() == out.shape(), 01647 "FFTWConvolvePlan::execute(): input and output must have the same shape."); 01648 01649 vigra_precondition(kernel.shape() == fourierArray.shape(), 01650 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan."); 01651 01652 Shape paddedShape = fftwCorrespondingShapeC2R(kernel.shape(), odd(in.shape(0))), 01653 diff = paddedShape - in.shape(), 01654 left = div(diff, MultiArrayIndex(2)), 01655 right = in.shape() + left; 01656 01657 vigra_precondition(paddedShape == realArray.shape(), 01658 "FFTWConvolvePlan::execute(): shape mismatch between input and plan."); 01659 01660 detail::fftEmbedArray(in, realArray); 01661 forward_plan.execute(realArray, fourierArray); 01662 01663 fourierKernel = kernel; 01664 moveDCToHalfspaceUpperLeft(fourierKernel); 01665 01666 fourierArray *= fourierKernel; 01667 01668 backward_plan.execute(fourierArray, realArray); 01669 01670 out = realArray.subarray(left, right); 01671 } 01672 01673 template <unsigned int N, class Real> 01674 template <class C1, class C2, class C3> 01675 void 01676 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01677 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 01678 MultiArrayView<N, FFTWComplex<Real>, C3> out) 01679 { 01680 vigra_precondition(in.shape() == out.shape(), 01681 "FFTWConvolvePlan::execute(): input and output must have the same shape."); 01682 01683 Shape paddedShape = fourierArray.shape(), 01684 diff = paddedShape - in.shape(), 01685 left = div(diff, MultiArrayIndex(2)), 01686 right = in.shape() + left; 01687 01688 if(useFourierKernel) 01689 { 01690 vigra_precondition(kernel.shape() == fourierArray.shape(), 01691 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan."); 01692 01693 fourierKernel = kernel; 01694 moveDCToUpperLeft(fourierKernel); 01695 } 01696 else 01697 { 01698 detail::fftEmbedKernel(kernel, fourierKernel); 01699 forward_plan.execute(fourierKernel, fourierKernel); 01700 } 01701 01702 detail::fftEmbedArray(in, fourierArray); 01703 forward_plan.execute(fourierArray, fourierArray); 01704 01705 fourierArray *= fourierKernel; 01706 01707 backward_plan.execute(fourierArray, fourierArray); 01708 01709 out = fourierArray.subarray(left, right); 01710 } 01711 01712 template <unsigned int N, class Real> 01713 template <class C1, class KernelIterator, class OutIterator> 01714 void 01715 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in, 01716 KernelIterator kernels, KernelIterator kernelsEnd, 01717 OutIterator outs, VigraFalseType /*useFourierKernel*/) 01718 { 01719 vigra_precondition(!useFourierKernel, 01720 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel."); 01721 01722 Shape kernelMax = checkShapes(in.shape(), kernels, kernelsEnd, outs), 01723 paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernelMax - Shape(1)), 01724 diff = paddedShape - in.shape(), 01725 left = div(diff, MultiArrayIndex(2)), 01726 right = in.shape() + left; 01727 01728 vigra_precondition(paddedShape == realArray.shape(), 01729 "FFTWConvolvePlan::executeMany(): shape mismatch between input and plan."); 01730 01731 detail::fftEmbedArray(in, realArray); 01732 forward_plan.execute(realArray, fourierArray); 01733 01734 for(; kernels != kernelsEnd; ++kernels, ++outs) 01735 { 01736 detail::fftEmbedKernel(*kernels, realKernel); 01737 forward_plan.execute(realKernel, fourierKernel); 01738 01739 fourierKernel *= fourierArray; 01740 01741 backward_plan.execute(fourierKernel, realKernel); 01742 01743 *outs = realKernel.subarray(left, right); 01744 } 01745 } 01746 01747 template <unsigned int N, class Real> 01748 template <class C1, class KernelIterator, class OutIterator> 01749 void 01750 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in, 01751 KernelIterator kernels, KernelIterator kernelsEnd, 01752 OutIterator outs, VigraTrueType /*useFourierKernel*/) 01753 { 01754 vigra_precondition(useFourierKernel, 01755 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel."); 01756 01757 Shape complexShape = checkShapesFourier(in.shape(), kernels, kernelsEnd, outs), 01758 paddedShape = fftwCorrespondingShapeC2R(complexShape, odd(in.shape(0))), 01759 diff = paddedShape - in.shape(), 01760 left = div(diff, MultiArrayIndex(2)), 01761 right = in.shape() + left; 01762 01763 vigra_precondition(complexShape == fourierArray.shape(), 01764 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between kernels and plan."); 01765 01766 vigra_precondition(paddedShape == realArray.shape(), 01767 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between input and plan."); 01768 01769 detail::fftEmbedArray(in, realArray); 01770 forward_plan.execute(realArray, fourierArray); 01771 01772 for(; kernels != kernelsEnd; ++kernels, ++outs) 01773 { 01774 fourierKernel = *kernels; 01775 moveDCToHalfspaceUpperLeft(fourierKernel); 01776 fourierKernel *= fourierArray; 01777 01778 backward_plan.execute(fourierKernel, realKernel); 01779 01780 *outs = realKernel.subarray(left, right); 01781 } 01782 } 01783 01784 template <unsigned int N, class Real> 01785 template <class C1, class KernelIterator, class OutIterator> 01786 void 01787 FFTWConvolvePlan<N, Real>::executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01788 KernelIterator kernels, KernelIterator kernelsEnd, 01789 OutIterator outs) 01790 { 01791 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray; 01792 typedef typename KernelArray::value_type KernelValue; 01793 typedef typename std::iterator_traits<OutIterator>::value_type OutArray; 01794 typedef typename OutArray::value_type OutValue; 01795 01796 vigra_precondition((IsSameType<KernelValue, Complex>::value), 01797 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type."); 01798 vigra_precondition((IsSameType<OutValue, Complex>::value), 01799 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type."); 01800 01801 Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs), 01802 diff = paddedShape - in.shape(), 01803 left = div(diff, MultiArrayIndex(2)), 01804 right = in.shape() + left; 01805 01806 detail::fftEmbedArray(in, fourierArray); 01807 forward_plan.execute(fourierArray, fourierArray); 01808 01809 for(; kernels != kernelsEnd; ++kernels, ++outs) 01810 { 01811 if(useFourierKernel) 01812 { 01813 fourierKernel = *kernels; 01814 moveDCToUpperLeft(fourierKernel); 01815 } 01816 else 01817 { 01818 detail::fftEmbedKernel(*kernels, fourierKernel); 01819 forward_plan.execute(fourierKernel, fourierKernel); 01820 } 01821 01822 fourierKernel *= fourierArray; 01823 01824 backward_plan.execute(fourierKernel, fourierKernel); 01825 01826 *outs = fourierKernel.subarray(left, right); 01827 } 01828 } 01829 01830 #endif // DOXYGEN 01831 01832 template <unsigned int N, class Real> 01833 template <class KernelIterator, class OutIterator> 01834 typename FFTWConvolvePlan<N, Real>::Shape 01835 FFTWConvolvePlan<N, Real>::checkShapes(Shape in, 01836 KernelIterator kernels, KernelIterator kernelsEnd, 01837 OutIterator outs) 01838 { 01839 vigra_precondition(kernels != kernelsEnd, 01840 "FFTWConvolvePlan::checkShapes(): empty kernel sequence."); 01841 01842 Shape kernelMax; 01843 for(; kernels != kernelsEnd; ++kernels, ++outs) 01844 { 01845 vigra_precondition(in == outs->shape(), 01846 "FFTWConvolvePlan::checkShapes(): shape mismatch between input and (one) output."); 01847 kernelMax = max(kernelMax, kernels->shape()); 01848 } 01849 vigra_precondition(prod(kernelMax) > 0, 01850 "FFTWConvolvePlan::checkShapes(): all kernels have size 0."); 01851 return kernelMax; 01852 } 01853 01854 template <unsigned int N, class Real> 01855 template <class KernelIterator, class OutIterator> 01856 typename FFTWConvolvePlan<N, Real>::Shape 01857 FFTWConvolvePlan<N, Real>::checkShapesFourier(Shape in, 01858 KernelIterator kernels, KernelIterator kernelsEnd, 01859 OutIterator outs) 01860 { 01861 vigra_precondition(kernels != kernelsEnd, 01862 "FFTWConvolvePlan::checkShapesFourier(): empty kernel sequence."); 01863 01864 Shape complexShape = kernels->shape(), 01865 paddedShape = fftwCorrespondingShapeC2R(complexShape); 01866 01867 for(unsigned int k=0; k<N; ++k) 01868 vigra_precondition(in[k] <= paddedShape[k], 01869 "FFTWConvolvePlan::checkShapesFourier(): kernels too small for given input."); 01870 01871 for(; kernels != kernelsEnd; ++kernels, ++outs) 01872 { 01873 vigra_precondition(in == outs->shape(), 01874 "FFTWConvolvePlan::checkShapesFourier(): shape mismatch between input and (one) output."); 01875 vigra_precondition(complexShape == kernels->shape(), 01876 "FFTWConvolvePlan::checkShapesFourier(): all kernels must have the same size."); 01877 } 01878 return complexShape; 01879 } 01880 01881 template <unsigned int N, class Real> 01882 template <class KernelIterator, class OutIterator> 01883 typename FFTWConvolvePlan<N, Real>::Shape 01884 FFTWConvolvePlan<N, Real>::checkShapesComplex(Shape in, 01885 KernelIterator kernels, KernelIterator kernelsEnd, 01886 OutIterator outs) 01887 { 01888 vigra_precondition(kernels != kernelsEnd, 01889 "FFTWConvolvePlan::checkShapesComplex(): empty kernel sequence."); 01890 01891 Shape kernelShape = kernels->shape(); 01892 for(; kernels != kernelsEnd; ++kernels, ++outs) 01893 { 01894 vigra_precondition(in == outs->shape(), 01895 "FFTWConvolvePlan::checkShapesComplex(): shape mismatch between input and (one) output."); 01896 if(useFourierKernel) 01897 { 01898 vigra_precondition(kernelShape == kernels->shape(), 01899 "FFTWConvolvePlan::checkShapesComplex(): Fourier domain kernels must have identical size."); 01900 } 01901 else 01902 { 01903 kernelShape = max(kernelShape, kernels->shape()); 01904 } 01905 } 01906 vigra_precondition(prod(kernelShape) > 0, 01907 "FFTWConvolvePlan::checkShapesComplex(): all kernels have size 0."); 01908 01909 if(useFourierKernel) 01910 { 01911 for(unsigned int k=0; k<N; ++k) 01912 vigra_precondition(in[k] <= kernelShape[k], 01913 "FFTWConvolvePlan::checkShapesComplex(): kernels too small for given input."); 01914 return kernelShape; 01915 } 01916 else 01917 { 01918 return fftwBestPaddedShape(in + kernelShape - Shape(1)); 01919 } 01920 } 01921 01922 01923 /********************************************************/ 01924 /* */ 01925 /* fourierTransform */ 01926 /* */ 01927 /********************************************************/ 01928 01929 template <unsigned int N, class Real, class C1, class C2> 01930 inline void 01931 fourierTransform(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01932 MultiArrayView<N, FFTWComplex<Real>, C2> out) 01933 { 01934 FFTWPlan<N, Real>(in, out, FFTW_FORWARD).execute(in, out); 01935 } 01936 01937 template <unsigned int N, class Real, class C1, class C2> 01938 inline void 01939 fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01940 MultiArrayView<N, FFTWComplex<Real>, C2> out) 01941 { 01942 FFTWPlan<N, Real>(in, out, FFTW_BACKWARD).execute(in, out); 01943 } 01944 01945 template <unsigned int N, class Real, class C1, class C2> 01946 void 01947 fourierTransform(MultiArrayView<N, Real, C1> in, 01948 MultiArrayView<N, FFTWComplex<Real>, C2> out) 01949 { 01950 if(in.shape() == out.shape()) 01951 { 01952 // copy the input array into the output and then perform an in-place FFT 01953 out = in; 01954 FFTWPlan<N, Real>(out, out, FFTW_FORWARD).execute(out, out); 01955 } 01956 else if(out.shape() == fftwCorrespondingShapeR2C(in.shape())) 01957 { 01958 FFTWPlan<N, Real>(in, out).execute(in, out); 01959 } 01960 else 01961 vigra_precondition(false, 01962 "fourierTransform(): shape mismatch between input and output."); 01963 } 01964 01965 template <unsigned int N, class Real, class C1, class C2> 01966 void 01967 fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01968 MultiArrayView<N, Real, C2> out) 01969 { 01970 vigra_precondition(in.shape() == fftwCorrespondingShapeR2C(out.shape()), 01971 "fourierTransformInverse(): shape mismatch between input and output."); 01972 FFTWPlan<N, Real>(in, out).execute(in, out); 01973 } 01974 01975 //@} 01976 01977 /** \addtogroup MultiArrayConvolutionFilters 01978 */ 01979 //@{ 01980 01981 /********************************************************/ 01982 /* */ 01983 /* convolveFFT */ 01984 /* */ 01985 /********************************************************/ 01986 01987 /** \brief Convolve an array with a kernel by means of the Fourier transform. 01988 01989 Thanks to the convolution theorem of Fourier theory, a convolution in the spatial domain 01990 is equivalent to a multiplication in the frequency domain. Thus, for certain kernels 01991 (especially large, non-separable ones), it is advantageous to perform the convolution by first 01992 transforming both array and kernel to the frequency domain, multiplying the frequency 01993 representations, and transforming the result back into the spatial domain. 01994 Some kernels have a much simpler definition in the frequency domain, so that they are readily 01995 computed there directly, avoiding Fourier transformation of those kernels. 01996 01997 The following functions implement various variants of FFT-based convolution: 01998 01999 <DL> 02000 <DT><b>convolveFFT</b><DD> Convolve a real-valued input array with a kernel such that the 02001 result is also real-valued. That is, the kernel is either provided 02002 as a real-valued array in the spatial domain, or as a 02003 complex-valued array in the Fourier domain, using the half-space format 02004 of the R2C Fourier transform (see below). 02005 <DT><b>convolveFFTMany</b><DD> Like <tt>convolveFFT</tt>, but you may provide many kernels at once 02006 (using an iterator pair specifying the kernel sequence). 02007 This has the advantage that the forward transform of the input array needs 02008 to be executed only once. 02009 <DT><b>convolveFFTComplex</b><DD> Convolve a complex-valued input array with a complex-valued kernel, 02010 resulting in a complex-valued output array. An additional flag is used to 02011 specify whether the kernel is defined in the spatial or frequency domain. 02012 <DT><b>convolveFFTComplexMany</b><DD> Like <tt>convolveFFTComplex</tt>, but you may provide many kernels at once 02013 (using an iterator pair specifying the kernel sequence). 02014 This has the advantage that the forward transform of the input array needs 02015 to be executed only once. 02016 </DL> 02017 02018 The output arrays must have the same shape as the input arrays. In the "Many" variants of the 02019 convolution functions, the kernels must all have the same shape. 02020 02021 The origin of the kernel is always assumed to be in the center of the kernel array (precisely, 02022 at the point <tt>floor(kernel.shape() / 2.0)</tt>, except when the half-space format is used, see below). 02023 The function \ref moveDCToUpperLeft() will be called internally to align the kernel with the transformed 02024 input as appropriate. 02025 02026 If a real input is combined with a real kernel, the kernel is automatically assumed to be defined 02027 in the spatial domain. If a real input is combined with a complex kernel, the kernel is assumed 02028 to be defined in the Fourier domain in half-space format. If the input array is complex, a flag 02029 <tt>fourierDomainKernel</tt> determines where the kernel is defined. 02030 02031 When the kernel is defined in the spatial domain, the convolution functions will automatically pad 02032 (enlarge) the input array by at least the kernel radius in each direction. The newly added space is 02033 filled according to reflective boundary conditions in order to minimize border artifacts during 02034 convolution. It is thus ensured that convolution in the Fourier domain yields the same results as 02035 convolution in the spatial domain (e.g. when \ref separableConvolveMultiArray() is called with the 02036 same kernel). A little further padding may be added to make sure that the padded array shape 02037 uses integers which have only small prime factors, because FFTW is then able to use the fastest 02038 possible algorithms. Any padding is automatically removed from the result arrays before the function 02039 returns. 02040 02041 When the kernel is defined in the frequency domain, it must be complex-valued, and its shape determines 02042 the shape of the Fourier representation (i.e. the input is padded according to the shape of the kernel). 02043 If we are going to perform a complex-valued convolution, the kernel must be defined for the entire 02044 frequency domain, and its shape directly determines the size of the FFT. 02045 02046 In contrast, a frequency domain kernel for a real-valued convolution must have symmetry properties 02047 that allow to drop half of the kernel coefficients, as in the 02048 <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C transform</a>. 02049 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 02050 logical shape of the frequency representation (and thus the size of the padded input). The origin 02051 of the kernel must be at the point 02052 <tt>(0, floor(fourier_shape[0] / 2.0), ..., floor(fourier_shape[N-1] / 2.0))</tt> 02053 (i.e. as in a regular kernel except for the first dimension). 02054 02055 The <tt>Real</tt> type in the declarations can be <tt>double</tt>, <tt>float</tt>, and 02056 <tt>long double</tt>. Your program must always link against <tt>libfftw3</tt>. If you use 02057 <tt>float</tt> or <tt>long double</tt> arrays, you must <i>additionally</i> link against 02058 <tt>libfftw3f</tt> and <tt>libfftw3l</tt> respectively. 02059 02060 The Fourier transform functions internally create <a href="http://www.fftw.org/doc/Using-Plans.html">FFTW plans</a> 02061 which control the algorithm details. The plans are creates with the flag <tt>FFTW_ESTIMATE</tt>, i.e. 02062 optimal settings are guessed or read from saved "wisdom" files. If you need more control over planning, 02063 you can use the class \ref FFTWConvolvePlan. 02064 02065 See also \ref applyFourierFilter() for corresponding functionality on the basis of the 02066 old image iterator interface. 02067 02068 <b> Declarations:</b> 02069 02070 Real-valued convolution with kernel in the spatial domain: 02071 \code 02072 namespace vigra { 02073 template <unsigned int N, class Real, class C1, class C2, class C3> 02074 void 02075 convolveFFT(MultiArrayView<N, Real, C1> in, 02076 MultiArrayView<N, Real, C2> kernel, 02077 MultiArrayView<N, Real, C3> out); 02078 } 02079 \endcode 02080 02081 Real-valued convolution with kernel in the Fourier domain (half-space format): 02082 \code 02083 namespace vigra { 02084 template <unsigned int N, class Real, class C1, class C2, class C3> 02085 void 02086 convolveFFT(MultiArrayView<N, Real, C1> in, 02087 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 02088 MultiArrayView<N, Real, C3> out); 02089 } 02090 \endcode 02091 02092 Series of real-valued convolutions with kernels in the spatial or Fourier domain 02093 (the kernel and out sequences must have the same length): 02094 \code 02095 namespace vigra { 02096 template <unsigned int N, class Real, class C1, 02097 class KernelIterator, class OutIterator> 02098 void 02099 convolveFFTMany(MultiArrayView<N, Real, C1> in, 02100 KernelIterator kernels, KernelIterator kernelsEnd, 02101 OutIterator outs); 02102 } 02103 \endcode 02104 02105 Complex-valued convolution (parameter <tt>fourierDomainKernel</tt> determines if 02106 the kernel is defined in the spatial or Fourier domain): 02107 \code 02108 namespace vigra { 02109 template <unsigned int N, class Real, class C1, class C2, class C3> 02110 void 02111 convolveFFTComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in, 02112 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 02113 MultiArrayView<N, FFTWComplex<Real>, C3> out, 02114 bool fourierDomainKernel); 02115 } 02116 \endcode 02117 02118 Series of complex-valued convolutions (parameter <tt>fourierDomainKernel</tt> 02119 determines if the kernels are defined in the spatial or Fourier domain, 02120 the kernel and out sequences must have the same length): 02121 \code 02122 namespace vigra { 02123 template <unsigned int N, class Real, class C1, 02124 class KernelIterator, class OutIterator> 02125 void 02126 convolveFFTComplexMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 02127 KernelIterator kernels, KernelIterator kernelsEnd, 02128 OutIterator outs, 02129 bool fourierDomainKernel); 02130 } 02131 \endcode 02132 02133 <b> Usage:</b> 02134 02135 <b>\#include</b> <vigra/multi_fft.hxx><br> 02136 Namespace: vigra 02137 02138 \code 02139 // convolve real array with a Gaussian (sigma=1) defined in the spatial domain 02140 // (implicitly uses padding by at least 4 pixels) 02141 MultiArray<2, double> src(Shape2(w, h)), dest(Shape2(w,h)); 02142 02143 MultiArray<2, double> spatial_kernel(Shape2(9, 9)); 02144 Gaussian<double> gauss(1.0); 02145 02146 for(int y=0; y<9; ++y) 02147 for(int x=0; x<9; ++x) 02148 spatial_kernel(x, y) = gauss(x-4.0)*gauss(y-4.0); 02149 02150 convolveFFT(src, spatial_kernel, dest); 02151 02152 // convolve real array with a Gaussian (sigma=1) defined in the Fourier domain 02153 // (uses no padding, because the kernel size corresponds to the input size) 02154 MultiArray<2, FFTWComplex<double> > fourier_kernel(fftwCorrespondingShapeR2C(src.shape())); 02155 int y0 = h / 2; 02156 02157 for(int y=0; y<fourier_kernel.shape(1); ++y) 02158 for(int x=0; x<fourier_kernel.shape(0); ++x) 02159 fourier_kernel(x, y) = exp(-0.5*sq(x / double(w))) * exp(-0.5*sq((y-y0)/double(h))); 02160 02161 convolveFFT(src, fourier_kernel, dest); 02162 \endcode 02163 */ 02164 doxygen_overloaded_function(template <...> void convolveFFT) 02165 02166 template <unsigned int N, class Real, class C1, class C2, class C3> 02167 void 02168 convolveFFT(MultiArrayView<N, Real, C1> in, 02169 MultiArrayView<N, Real, C2> kernel, 02170 MultiArrayView<N, Real, C3> out) 02171 { 02172 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out); 02173 } 02174 02175 template <unsigned int N, class Real, class C1, class C2, class C3> 02176 void 02177 convolveFFT(MultiArrayView<N, Real, C1> in, 02178 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 02179 MultiArrayView<N, Real, C3> out) 02180 { 02181 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out); 02182 } 02183 02184 /** \brief Convolve a complex-valued array by means of the Fourier transform. 02185 02186 See \ref convolveFFT() for details. 02187 */ 02188 doxygen_overloaded_function(template <...> void convolveFFTComplex) 02189 02190 template <unsigned int N, class Real, class C1, class C2, class C3> 02191 void 02192 convolveFFTComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in, 02193 MultiArrayView<N, FFTWComplex<Real>, C2> kernel, 02194 MultiArrayView<N, FFTWComplex<Real>, C3> out, 02195 bool fourierDomainKernel) 02196 { 02197 FFTWConvolvePlan<N, Real>(in, kernel, out, fourierDomainKernel).execute(in, kernel, out); 02198 } 02199 02200 /** \brief Convolve a real-valued array with a sequence of kernels by means of the Fourier transform. 02201 02202 See \ref convolveFFT() for details. 02203 */ 02204 doxygen_overloaded_function(template <...> void convolveFFTMany) 02205 02206 template <unsigned int N, class Real, class C1, 02207 class KernelIterator, class OutIterator> 02208 void 02209 convolveFFTMany(MultiArrayView<N, Real, C1> in, 02210 KernelIterator kernels, KernelIterator kernelsEnd, 02211 OutIterator outs) 02212 { 02213 FFTWConvolvePlan<N, Real> plan; 02214 plan.initMany(in, kernels, kernelsEnd, outs); 02215 plan.executeMany(in, kernels, kernelsEnd, outs); 02216 } 02217 02218 /** \brief Convolve a complex-valued array with a sequence of kernels by means of the Fourier transform. 02219 02220 See \ref convolveFFT() for details. 02221 */ 02222 doxygen_overloaded_function(template <...> void convolveFFTComplexMany) 02223 02224 template <unsigned int N, class Real, class C1, 02225 class KernelIterator, class OutIterator> 02226 void 02227 convolveFFTComplexMany(MultiArrayView<N, FFTWComplex<Real>, C1> in, 02228 KernelIterator kernels, KernelIterator kernelsEnd, 02229 OutIterator outs, 02230 bool fourierDomainKernel) 02231 { 02232 FFTWConvolvePlan<N, Real> plan; 02233 plan.initMany(in, kernels, kernelsEnd, outs, fourierDomainKernel); 02234 plan.executeMany(in, kernels, kernelsEnd, outs); 02235 } 02236 02237 //@} 02238 02239 } // namespace vigra 02240 02241 #endif // VIGRA_MULTI_FFT_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|