[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/multi_fft.hxx VIGRA

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)