vigra - vigra/multi_fft.hxx Source File

[ 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) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)