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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (20 Sep 2011)