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

multi_fft.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2009-2010 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_MULTI_FFT_HXX
37 #define VIGRA_MULTI_FFT_HXX
38 
39 #include "fftw3.hxx"
40 #include "multi_array.hxx"
41 #include "navigator.hxx"
42 #include "copyimage.hxx"
43 
44 namespace vigra {
45 
46 /********************************************************/
47 /* */
48 /* Fourier Transform */
49 /* */
50 /********************************************************/
51 
52 /** \addtogroup FourierTransform
53 */
54 //@{
55 
56 namespace detail {
57 
58 template <unsigned int N, class T, class C>
59 void moveDCToCenterImpl(MultiArrayView<N, T, C> a, unsigned int startDimension)
60 {
61  typedef typename MultiArrayView<N, T, C>::traverser Traverser;
62  typedef MultiArrayNavigator<Traverser, N> Navigator;
63  typedef typename Navigator::iterator Iterator;
64 
65  for(unsigned int d = startDimension; d < N; ++d)
66  {
67  Navigator nav(a.traverser_begin(), a.shape(), d);
68 
69  for( ; nav.hasMore(); nav++ )
70  {
71  Iterator i = nav.begin();
72  int s = nav.end() - i;
73  int s2 = s/2;
74 
75  if(even(s))
76  {
77  for(int k=0; k<s2; ++k)
78  {
79  std::swap(i[k], i[k+s2]);
80  }
81  }
82  else
83  {
84  T v = i[0];
85  for(int k=0; k<s2; ++k)
86  {
87  i[k] = i[k+s2+1];
88  i[k+s2+1] = i[k+1];
89  }
90  i[s2] = v;
91  }
92  }
93  }
94 }
95 
96 template <unsigned int N, class T, class C>
97 void moveDCToUpperLeftImpl(MultiArrayView<N, T, C> a, unsigned int startDimension)
98 {
99  typedef typename MultiArrayView<N, T, C>::traverser Traverser;
100  typedef MultiArrayNavigator<Traverser, N> Navigator;
101  typedef typename Navigator::iterator Iterator;
102 
103  for(unsigned int d = startDimension; d < N; ++d)
104  {
105  Navigator nav(a.traverser_begin(), a.shape(), d);
106 
107  for( ; nav.hasMore(); nav++ )
108  {
109  Iterator i = nav.begin();
110  int s = nav.end() - i;
111  int s2 = s/2;
112 
113  if(even(s))
114  {
115  for(int k=0; k<s2; ++k)
116  {
117  std::swap(i[k], i[k+s2]);
118  }
119  }
120  else
121  {
122  T v = i[s2];
123  for(int k=s2; k>0; --k)
124  {
125  i[k] = i[k+s2];
126  i[k+s2] = i[k-1];
127  }
128  i[0] = v;
129  }
130  }
131  }
132 }
133 
134 } // namespace detail
135 
136 /********************************************************/
137 /* */
138 /* moveDCToCenter */
139 /* */
140 /********************************************************/
141 
142 template <unsigned int N, class T, class C>
143 inline void moveDCToCenter(MultiArrayView<N, T, C> a)
144 {
145  detail::moveDCToCenterImpl(a, 0);
146 }
147 
148 template <unsigned int N, class T, class C>
149 inline void moveDCToUpperLeft(MultiArrayView<N, T, C> a)
150 {
151  detail::moveDCToUpperLeftImpl(a, 0);
152 }
153 
154 template <unsigned int N, class T, class C>
155 inline void moveDCToHalfspaceCenter(MultiArrayView<N, T, C> a)
156 {
157  detail::moveDCToCenterImpl(a, 1);
158 }
159 
160 template <unsigned int N, class T, class C>
161 inline void moveDCToHalfspaceUpperLeft(MultiArrayView<N, T, C> a)
162 {
163  detail::moveDCToUpperLeftImpl(a, 1);
164 }
165 
166 namespace detail
167 {
168 
169 inline fftw_plan
170 fftwPlanCreate(unsigned int N, int* shape,
171  FFTWComplex<double> * in, int* instrides, int instep,
172  FFTWComplex<double> * out, int* outstrides, int outstep,
173  int sign, unsigned int planner_flags)
174 {
175  return fftw_plan_many_dft(N, shape, 1,
176  (fftw_complex *)in, instrides, instep, 0,
177  (fftw_complex *)out, outstrides, outstep, 0,
178  sign, planner_flags);
179 }
180 
181 inline fftw_plan
182 fftwPlanCreate(unsigned int N, int* shape,
183  double * in, int* instrides, int instep,
184  FFTWComplex<double> * out, int* outstrides, int outstep,
185  int /*sign is ignored*/, unsigned int planner_flags)
186 {
187  return fftw_plan_many_dft_r2c(N, shape, 1,
188  in, instrides, instep, 0,
189  (fftw_complex *)out, outstrides, outstep, 0,
190  planner_flags);
191 }
192 
193 inline fftw_plan
194 fftwPlanCreate(unsigned int N, int* shape,
195  FFTWComplex<double> * in, int* instrides, int instep,
196  double * out, int* outstrides, int outstep,
197  int /*sign is ignored*/, unsigned int planner_flags)
198 {
199  return fftw_plan_many_dft_c2r(N, shape, 1,
200  (fftw_complex *)in, instrides, instep, 0,
201  out, outstrides, outstep, 0,
202  planner_flags);
203 }
204 
205 inline fftwf_plan
206 fftwPlanCreate(unsigned int N, int* shape,
207  FFTWComplex<float> * in, int* instrides, int instep,
208  FFTWComplex<float> * out, int* outstrides, int outstep,
209  int sign, unsigned int planner_flags)
210 {
211  return fftwf_plan_many_dft(N, shape, 1,
212  (fftwf_complex *)in, instrides, instep, 0,
213  (fftwf_complex *)out, outstrides, outstep, 0,
214  sign, planner_flags);
215 }
216 
217 inline fftwf_plan
218 fftwPlanCreate(unsigned int N, int* shape,
219  float * in, int* instrides, int instep,
220  FFTWComplex<float> * out, int* outstrides, int outstep,
221  int /*sign is ignored*/, unsigned int planner_flags)
222 {
223  return fftwf_plan_many_dft_r2c(N, shape, 1,
224  in, instrides, instep, 0,
225  (fftwf_complex *)out, outstrides, outstep, 0,
226  planner_flags);
227 }
228 
229 inline fftwf_plan
230 fftwPlanCreate(unsigned int N, int* shape,
231  FFTWComplex<float> * in, int* instrides, int instep,
232  float * out, int* outstrides, int outstep,
233  int /*sign is ignored*/, unsigned int planner_flags)
234 {
235  return fftwf_plan_many_dft_c2r(N, shape, 1,
236  (fftwf_complex *)in, instrides, instep, 0,
237  out, outstrides, outstep, 0,
238  planner_flags);
239 }
240 
241 inline fftwl_plan
242 fftwPlanCreate(unsigned int N, int* shape,
243  FFTWComplex<long double> * in, int* instrides, int instep,
244  FFTWComplex<long double> * out, int* outstrides, int outstep,
245  int sign, unsigned int planner_flags)
246 {
247  return fftwl_plan_many_dft(N, shape, 1,
248  (fftwl_complex *)in, instrides, instep, 0,
249  (fftwl_complex *)out, outstrides, outstep, 0,
250  sign, planner_flags);
251 }
252 
253 inline fftwl_plan
254 fftwPlanCreate(unsigned int N, int* shape,
255  long double * in, int* instrides, int instep,
256  FFTWComplex<long double> * out, int* outstrides, int outstep,
257  int /*sign is ignored*/, unsigned int planner_flags)
258 {
259  return fftwl_plan_many_dft_r2c(N, shape, 1,
260  in, instrides, instep, 0,
261  (fftwl_complex *)out, outstrides, outstep, 0,
262  planner_flags);
263 }
264 
265 inline fftwl_plan
266 fftwPlanCreate(unsigned int N, int* shape,
267  FFTWComplex<long double> * in, int* instrides, int instep,
268  long double * out, int* outstrides, int outstep,
269  int /*sign is ignored*/, unsigned int planner_flags)
270 {
271  return fftwl_plan_many_dft_c2r(N, shape, 1,
272  (fftwl_complex *)in, instrides, instep, 0,
273  out, outstrides, outstep, 0,
274  planner_flags);
275 }
276 
277 inline void fftwPlanDestroy(fftw_plan plan)
278 {
279  if(plan != 0)
280  fftw_destroy_plan(plan);
281 }
282 
283 inline void fftwPlanDestroy(fftwf_plan plan)
284 {
285  if(plan != 0)
286  fftwf_destroy_plan(plan);
287 }
288 
289 inline void fftwPlanDestroy(fftwl_plan plan)
290 {
291  if(plan != 0)
292  fftwl_destroy_plan(plan);
293 }
294 
295 inline void
296 fftwPlanExecute(fftw_plan plan)
297 {
298  fftw_execute(plan);
299 }
300 
301 inline void
302 fftwPlanExecute(fftwf_plan plan)
303 {
304  fftwf_execute(plan);
305 }
306 
307 inline void
308 fftwPlanExecute(fftwl_plan plan)
309 {
310  fftwl_execute(plan);
311 }
312 
313 inline void
314 fftwPlanExecute(fftw_plan plan, FFTWComplex<double> * in, FFTWComplex<double> * out)
315 {
316  fftw_execute_dft(plan, (fftw_complex *)in, (fftw_complex *)out);
317 }
318 
319 inline void
320 fftwPlanExecute(fftw_plan plan, double * in, FFTWComplex<double> * out)
321 {
322  fftw_execute_dft_r2c(plan, in, (fftw_complex *)out);
323 }
324 
325 inline void
326 fftwPlanExecute(fftw_plan plan, FFTWComplex<double> * in, double * out)
327 {
328  fftw_execute_dft_c2r(plan, (fftw_complex *)in, out);
329 }
330 
331 inline void
332 fftwPlanExecute(fftwf_plan plan, FFTWComplex<float> * in, FFTWComplex<float> * out)
333 {
334  fftwf_execute_dft(plan, (fftwf_complex *)in, (fftwf_complex *)out);
335 }
336 
337 inline void
338 fftwPlanExecute(fftwf_plan plan, float * in, FFTWComplex<float> * out)
339 {
340  fftwf_execute_dft_r2c(plan, in, (fftwf_complex *)out);
341 }
342 
343 inline void
344 fftwPlanExecute(fftwf_plan plan, FFTWComplex<float> * in, float * out)
345 {
346  fftwf_execute_dft_c2r(plan, (fftwf_complex *)in, out);
347 }
348 
349 inline void
350 fftwPlanExecute(fftwl_plan plan, FFTWComplex<long double> * in, FFTWComplex<long double> * out)
351 {
352  fftwl_execute_dft(plan, (fftwl_complex *)in, (fftwl_complex *)out);
353 }
354 
355 inline void
356 fftwPlanExecute(fftwl_plan plan, long double * in, FFTWComplex<long double> * out)
357 {
358  fftwl_execute_dft_r2c(plan, in, (fftwl_complex *)out);
359 }
360 
361 inline void
362 fftwPlanExecute(fftwl_plan plan, FFTWComplex<long double> * in, long double * out)
363 {
364  fftwl_execute_dft_c2r(plan, (fftwl_complex *)in, out);
365 }
366 
367 inline
368 int fftwPaddingSize(int s)
369 {
370  // Image sizes where FFTW is fast. The list contains all
371  // numbers less than 100000 whose prime decomposition is of the form
372  // 2^a*3^b*5^c*7^d*11^e*13^f, where e+f is either 0 or 1, and the
373  // other exponents are arbitrary
374  static const int size = 1330;
375  static int goodSizes[size] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
376  18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 35, 36, 39, 40, 42, 44, 45, 48,
377  49, 50, 52, 54, 55, 56, 60, 63, 64, 65, 66, 70, 72, 75, 77, 78, 80, 81,
378  84, 88, 90, 91, 96, 98, 99, 100, 104, 105, 108, 110, 112, 117, 120, 125,
379  126, 128, 130, 132, 135, 140, 144, 147, 150, 154, 156, 160, 162, 165,
380  168, 175, 176, 180, 182, 189, 192, 195, 196, 198, 200, 208, 210, 216,
381  220, 224, 225, 231, 234, 240, 243, 245, 250, 252, 256, 260, 264, 270,
382  273, 275, 280, 288, 294, 297, 300, 308, 312, 315, 320, 324, 325, 330,
383  336, 343, 350, 351, 352, 360, 364, 375, 378, 384, 385, 390, 392, 396,
384  400, 405, 416, 420, 432, 440, 441, 448, 450, 455, 462, 468, 480, 486,
385  490, 495, 500, 504, 512, 520, 525, 528, 539, 540, 546, 550, 560, 567,
386  576, 585, 588, 594, 600, 616, 624, 625, 630, 637, 640, 648, 650, 660,
387  672, 675, 686, 693, 700, 702, 704, 720, 728, 729, 735, 750, 756, 768,
388  770, 780, 784, 792, 800, 810, 819, 825, 832, 840, 864, 875, 880, 882,
389  891, 896, 900, 910, 924, 936, 945, 960, 972, 975, 980, 990, 1000, 1008,
390  1024, 1029, 1040, 1050, 1053, 1056, 1078, 1080, 1092, 1100, 1120, 1125,
391  1134, 1152, 1155, 1170, 1176, 1188, 1200, 1215, 1225, 1232, 1248, 1250,
392  1260, 1274, 1280, 1296, 1300, 1320, 1323, 1344, 1350, 1365, 1372, 1375,
393  1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470, 1485, 1500, 1512, 1536,
394  1540, 1560, 1568, 1575, 1584, 1600, 1617, 1620, 1625, 1638, 1650, 1664,
395  1680, 1701, 1715, 1728, 1750, 1755, 1760, 1764, 1782, 1792, 1800, 1820,
396  1848, 1872, 1875, 1890, 1911, 1920, 1925, 1944, 1950, 1960, 1980, 2000,
397  2016, 2025, 2048, 2058, 2079, 2080, 2100, 2106, 2112, 2156, 2160, 2184,
398  2187, 2200, 2205, 2240, 2250, 2268, 2275, 2304, 2310, 2340, 2352, 2376,
399  2400, 2401, 2430, 2450, 2457, 2464, 2475, 2496, 2500, 2520, 2548, 2560,
400  2592, 2600, 2625, 2640, 2646, 2673, 2688, 2695, 2700, 2730, 2744, 2750,
401  2772, 2800, 2808, 2816, 2835, 2880, 2912, 2916, 2925, 2940, 2970, 3000,
402  3024, 3072, 3080, 3087, 3120, 3125, 3136, 3150, 3159, 3168, 3185, 3200,
403  3234, 3240, 3250, 3276, 3300, 3328, 3360, 3375, 3402, 3430, 3456, 3465,
404  3500, 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3645, 3675, 3696, 3744,
405  3750, 3773, 3780, 3822, 3840, 3850, 3888, 3900, 3920, 3960, 3969, 4000,
406  4032, 4050, 4095, 4096, 4116, 4125, 4158, 4160, 4200, 4212, 4224, 4312,
407  4320, 4368, 4374, 4375, 4400, 4410, 4455, 4459, 4480, 4500, 4536, 4550,
408  4608, 4620, 4680, 4704, 4725, 4752, 4800, 4802, 4851, 4860, 4875, 4900,
409  4914, 4928, 4950, 4992, 5000, 5040, 5096, 5103, 5120, 5145, 5184, 5200,
410  5250, 5265, 5280, 5292, 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544,
411  5600, 5616, 5625, 5632, 5670, 5733, 5760, 5775, 5824, 5832, 5850, 5880,
412  5940, 6000, 6048, 6075, 6125, 6144, 6160, 6174, 6237, 6240, 6250, 6272,
413  6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6561, 6600, 6615,
414  6656, 6720, 6750, 6804, 6825, 6860, 6875, 6912, 6930, 7000, 7020, 7040,
415  7056, 7128, 7168, 7200, 7203, 7280, 7290, 7350, 7371, 7392, 7425, 7488,
416  7500, 7546, 7560, 7644, 7680, 7700, 7776, 7800, 7840, 7875, 7920, 7938,
417  8000, 8019, 8064, 8085, 8100, 8125, 8190, 8192, 8232, 8250, 8316, 8320,
418  8400, 8424, 8448, 8505, 8575, 8624, 8640, 8736, 8748, 8750, 8775, 8800,
419  8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9261, 9360, 9375,
420  9408, 9450, 9477, 9504, 9555, 9600, 9604, 9625, 9702, 9720, 9750, 9800,
421  9828, 9856, 9900, 9984, 10000, 10080, 10125, 10192, 10206, 10240, 10290,
422  10368, 10395, 10400, 10500, 10530, 10560, 10584, 10692, 10752, 10780,
423  10800, 10920, 10935, 10976, 11000, 11025, 11088, 11200, 11232, 11250,
424  11264, 11319, 11340, 11375, 11466, 11520, 11550, 11648, 11664, 11700,
425  11760, 11880, 11907, 12000, 12005, 12096, 12150, 12250, 12285, 12288,
426  12320, 12348, 12375, 12474, 12480, 12500, 12544, 12600, 12636, 12672,
427  12740, 12800, 12936, 12960, 13000, 13104, 13122, 13125, 13200, 13230,
428  13312, 13365, 13377, 13440, 13475, 13500, 13608, 13650, 13720, 13750,
429  13824, 13860, 14000, 14040, 14080, 14112, 14175, 14256, 14336, 14400,
430  14406, 14553, 14560, 14580, 14625, 14700, 14742, 14784, 14850, 14976,
431  15000, 15092, 15120, 15288, 15309, 15360, 15400, 15435, 15552, 15600,
432  15625, 15680, 15750, 15795, 15840, 15876, 15925, 16000, 16038, 16128,
433  16170, 16200, 16250, 16380, 16384, 16464, 16500, 16632, 16640, 16800,
434  16807, 16848, 16875, 16896, 17010, 17150, 17199, 17248, 17280, 17325,
435  17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000,
436  18144, 18200, 18225, 18375, 18432, 18480, 18522, 18711, 18720, 18750,
437  18816, 18865, 18900, 18954, 19008, 19110, 19200, 19208, 19250, 19404,
438  19440, 19500, 19600, 19656, 19683, 19712, 19800, 19845, 19968, 20000,
439  20160, 20250, 20384, 20412, 20475, 20480, 20580, 20625, 20736, 20790,
440  20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560, 21600, 21609,
441  21840, 21870, 21875, 21952, 22000, 22050, 22113, 22176, 22275, 22295,
442  22400, 22464, 22500, 22528, 22638, 22680, 22750, 22932, 23040, 23100,
443  23296, 23328, 23400, 23520, 23625, 23760, 23814, 24000, 24010, 24057,
444  24192, 24255, 24300, 24375, 24500, 24570, 24576, 24640, 24696, 24750,
445  24948, 24960, 25000, 25088, 25200, 25272, 25344, 25480, 25515, 25600,
446  25725, 25872, 25920, 26000, 26208, 26244, 26250, 26325, 26400, 26411,
447  26460, 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440,
448  27500, 27648, 27720, 27783, 28000, 28080, 28125, 28160, 28224, 28350,
449  28431, 28512, 28665, 28672, 28800, 28812, 28875, 29106, 29120, 29160,
450  29250, 29400, 29484, 29568, 29700, 29952, 30000, 30184, 30240, 30375,
451  30576, 30618, 30625, 30720, 30800, 30870, 31104, 31185, 31200, 31213,
452  31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256,
453  32340, 32400, 32500, 32760, 32768, 32805, 32928, 33000, 33075, 33264,
454  33280, 33600, 33614, 33696, 33750, 33792, 33957, 34020, 34125, 34300,
455  34375, 34398, 34496, 34560, 34650, 34944, 34992, 35000, 35100, 35200,
456  35280, 35640, 35672, 35721, 35840, 36000, 36015, 36288, 36400, 36450,
457  36750, 36855, 36864, 36960, 37044, 37125, 37422, 37440, 37500, 37632,
458  37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808, 38880,
459  39000, 39200, 39312, 39366, 39375, 39424, 39600, 39690, 39936, 40000,
460  40095, 40131, 40320, 40425, 40500, 40625, 40768, 40824, 40950, 40960,
461  41160, 41250, 41472, 41580, 41600, 42000, 42120, 42240, 42336, 42525,
462  42768, 42875, 43008, 43120, 43200, 43218, 43659, 43680, 43740, 43750,
463  43875, 43904, 44000, 44100, 44226, 44352, 44550, 44590, 44800, 44928,
464  45000, 45056, 45276, 45360, 45500, 45864, 45927, 46080, 46200, 46305,
465  46592, 46656, 46800, 46875, 47040, 47250, 47385, 47520, 47628, 47775,
466  48000, 48020, 48114, 48125, 48384, 48510, 48600, 48750, 49000, 49140,
467  49152, 49280, 49392, 49500, 49896, 49920, 50000, 50176, 50400, 50421,
468  50544, 50625, 50688, 50960, 51030, 51200, 51450, 51597, 51744, 51840,
469  51975, 52000, 52416, 52488, 52500, 52650, 52800, 52822, 52920, 53248,
470  53460, 53508, 53760, 53900, 54000, 54432, 54600, 54675, 54880, 55000,
471  55125, 55296, 55440, 55566, 56000, 56133, 56160, 56250, 56320, 56448,
472  56595, 56700, 56862, 56875, 57024, 57330, 57344, 57600, 57624, 57750,
473  58212, 58240, 58320, 58500, 58800, 58968, 59049, 59136, 59400, 59535,
474  59904, 60000, 60025, 60368, 60480, 60750, 61152, 61236, 61250, 61425,
475  61440, 61600, 61740, 61875, 62208, 62370, 62400, 62426, 62500, 62720,
476  63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512, 64680, 64800,
477  64827, 65000, 65520, 65536, 65610, 65625, 65856, 66000, 66150, 66339,
478  66528, 66560, 66825, 66885, 67200, 67228, 67375, 67392, 67500, 67584,
479  67914, 68040, 68250, 68600, 68750, 68796, 68992, 69120, 69300, 69888,
480  69984, 70000, 70200, 70400, 70560, 70875, 71280, 71344, 71442, 71680,
481  72000, 72030, 72171, 72576, 72765, 72800, 72900, 73125, 73500, 73710,
482  73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264, 75460, 75600,
483  75816, 76032, 76440, 76545, 76800, 76832, 77000, 77175, 77616, 77760,
484  78000, 78125, 78400, 78624, 78732, 78750, 78848, 78975, 79200, 79233,
485  79380, 79625, 79872, 80000, 80190, 80262, 80640, 80850, 81000, 81250,
486  81536, 81648, 81900, 81920, 82320, 82500, 82944, 83160, 83200, 83349,
487  84000, 84035, 84240, 84375, 84480, 84672, 85050, 85293, 85536, 85750,
488  85995, 86016, 86240, 86400, 86436, 86625, 87318, 87360, 87480, 87500,
489  87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856,
490  90000, 90112, 90552, 90720, 91000, 91125, 91728, 91854, 91875, 92160,
491  92400, 92610, 93184, 93312, 93555, 93600, 93639, 93750, 94080, 94325,
492  94500, 94770, 95040, 95256, 95550, 96000, 96040, 96228, 96250, 96768,
493  97020, 97200, 97500, 98000, 98280, 98304, 98415, 98560, 98784, 99000,
494  99225, 99792, 99840 };
495 
496  if(s <= 0 || s > goodSizes[size-1])
497  return s;
498  return *std::upper_bound(goodSizes, goodSizes+size, s, std::less_equal<int>());
499 }
500 
501 inline
502 int fftwEvenPaddingSize(int s)
503 {
504  // Image sizes where FFTW is fast. The list contains all even
505  // numbers less than 100000 whose prime decomposition is of the form
506  // 2^a*3^b*5^c*7^d*11^e*13^f, where a >= 1, e+f is either 0 or 1, and the
507  // other exponents are arbitrary
508  static const int size = 1063;
509  static int goodSizes[size] = { 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22,
510  24, 26, 28, 30, 32, 36, 40, 42, 44, 48, 50, 52, 54, 56, 60, 64, 66, 70,
511  72, 78, 80, 84, 88, 90, 96, 98, 100, 104, 108, 110, 112, 120, 126, 128,
512  130, 132, 140, 144, 150, 154, 156, 160, 162, 168, 176, 180, 182, 192,
513  196, 198, 200, 208, 210, 216, 220, 224, 234, 240, 250, 252, 256, 260,
514  264, 270, 280, 288, 294, 300, 308, 312, 320, 324, 330, 336, 350, 352,
515  360, 364, 378, 384, 390, 392, 396, 400, 416, 420, 432, 440, 448, 450,
516  462, 468, 480, 486, 490, 500, 504, 512, 520, 528, 540, 546, 550, 560,
517  576, 588, 594, 600, 616, 624, 630, 640, 648, 650, 660, 672, 686, 700,
518  702, 704, 720, 728, 750, 756, 768, 770, 780, 784, 792, 800, 810, 832,
519  840, 864, 880, 882, 896, 900, 910, 924, 936, 960, 972, 980, 990, 1000,
520  1008, 1024, 1040, 1050, 1056, 1078, 1080, 1092, 1100, 1120, 1134, 1152,
521  1170, 1176, 1188, 1200, 1232, 1248, 1250, 1260, 1274, 1280, 1296, 1300,
522  1320, 1344, 1350, 1372, 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470,
523  1500, 1512, 1536, 1540, 1560, 1568, 1584, 1600, 1620, 1638, 1650, 1664,
524  1680, 1728, 1750, 1760, 1764, 1782, 1792, 1800, 1820, 1848, 1872, 1890,
525  1920, 1944, 1950, 1960, 1980, 2000, 2016, 2048, 2058, 2080, 2100, 2106,
526  2112, 2156, 2160, 2184, 2200, 2240, 2250, 2268, 2304, 2310, 2340, 2352,
527  2376, 2400, 2430, 2450, 2464, 2496, 2500, 2520, 2548, 2560, 2592, 2600,
528  2640, 2646, 2688, 2700, 2730, 2744, 2750, 2772, 2800, 2808, 2816, 2880,
529  2912, 2916, 2940, 2970, 3000, 3024, 3072, 3080, 3120, 3136, 3150, 3168,
530  3200, 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3402, 3430, 3456, 3500,
531  3510, 3520, 3528, 3564, 3584, 3600, 3640, 3696, 3744, 3750, 3780, 3822,
532  3840, 3850, 3888, 3900, 3920, 3960, 4000, 4032, 4050, 4096, 4116, 4158,
533  4160, 4200, 4212, 4224, 4312, 4320, 4368, 4374, 4400, 4410, 4480, 4500,
534  4536, 4550, 4608, 4620, 4680, 4704, 4752, 4800, 4802, 4860, 4900, 4914,
535  4928, 4950, 4992, 5000, 5040, 5096, 5120, 5184, 5200, 5250, 5280, 5292,
536  5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544, 5600, 5616, 5632, 5670,
537  5760, 5824, 5832, 5850, 5880, 5940, 6000, 6048, 6144, 6160, 6174, 6240,
538  6250, 6272, 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6600,
539  6656, 6720, 6750, 6804, 6860, 6912, 6930, 7000, 7020, 7040, 7056, 7128,
540  7168, 7200, 7280, 7290, 7350, 7392, 7488, 7500, 7546, 7560, 7644, 7680,
541  7700, 7776, 7800, 7840, 7920, 7938, 8000, 8064, 8100, 8190, 8192, 8232,
542  8250, 8316, 8320, 8400, 8424, 8448, 8624, 8640, 8736, 8748, 8750, 8800,
543  8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9360, 9408, 9450,
544  9504, 9600, 9604, 9702, 9720, 9750, 9800, 9828, 9856, 9900, 9984, 10000,
545  10080, 10192, 10206, 10240, 10290, 10368, 10400, 10500, 10530, 10560,
546  10584, 10692, 10752, 10780, 10800, 10920, 10976, 11000, 11088, 11200,
547  11232, 11250, 11264, 11340, 11466, 11520, 11550, 11648, 11664, 11700,
548  11760, 11880, 12000, 12096, 12150, 12250, 12288, 12320, 12348, 12474,
549  12480, 12500, 12544, 12600, 12636, 12672, 12740, 12800, 12936, 12960,
550  13000, 13104, 13122, 13200, 13230, 13312, 13440, 13500, 13608, 13650,
551  13720, 13750, 13824, 13860, 14000, 14040, 14080, 14112, 14256, 14336,
552  14400, 14406, 14560, 14580, 14700, 14742, 14784, 14850, 14976, 15000,
553  15092, 15120, 15288, 15360, 15400, 15552, 15600, 15680, 15750, 15840,
554  15876, 16000, 16038, 16128, 16170, 16200, 16250, 16380, 16384, 16464,
555  16500, 16632, 16640, 16800, 16848, 16896, 17010, 17150, 17248, 17280,
556  17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000,
557  18144, 18200, 18432, 18480, 18522, 18720, 18750, 18816, 18900, 18954,
558  19008, 19110, 19200, 19208, 19250, 19404, 19440, 19500, 19600, 19656,
559  19712, 19800, 19968, 20000, 20160, 20250, 20384, 20412, 20480, 20580,
560  20736, 20790, 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560,
561  21600, 21840, 21870, 21952, 22000, 22050, 22176, 22400, 22464, 22500,
562  22528, 22638, 22680, 22750, 22932, 23040, 23100, 23296, 23328, 23400,
563  23520, 23760, 23814, 24000, 24010, 24192, 24300, 24500, 24570, 24576,
564  24640, 24696, 24750, 24948, 24960, 25000, 25088, 25200, 25272, 25344,
565  25480, 25600, 25872, 25920, 26000, 26208, 26244, 26250, 26400, 26460,
566  26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440, 27500,
567  27648, 27720, 28000, 28080, 28160, 28224, 28350, 28512, 28672, 28800,
568  28812, 29106, 29120, 29160, 29250, 29400, 29484, 29568, 29700, 29952,
569  30000, 30184, 30240, 30576, 30618, 30720, 30800, 30870, 31104, 31200,
570  31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256,
571  32340, 32400, 32500, 32760, 32768, 32928, 33000, 33264, 33280, 33600,
572  33614, 33696, 33750, 33792, 34020, 34300, 34398, 34496, 34560, 34650,
573  34944, 34992, 35000, 35100, 35200, 35280, 35640, 35672, 35840, 36000,
574  36288, 36400, 36450, 36750, 36864, 36960, 37044, 37422, 37440, 37500,
575  37632, 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808,
576  38880, 39000, 39200, 39312, 39366, 39424, 39600, 39690, 39936, 40000,
577  40320, 40500, 40768, 40824, 40950, 40960, 41160, 41250, 41472, 41580,
578  41600, 42000, 42120, 42240, 42336, 42768, 43008, 43120, 43200, 43218,
579  43680, 43740, 43750, 43904, 44000, 44100, 44226, 44352, 44550, 44590,
580  44800, 44928, 45000, 45056, 45276, 45360, 45500, 45864, 46080, 46200,
581  46592, 46656, 46800, 47040, 47250, 47520, 47628, 48000, 48020, 48114,
582  48384, 48510, 48600, 48750, 49000, 49140, 49152, 49280, 49392, 49500,
583  49896, 49920, 50000, 50176, 50400, 50544, 50688, 50960, 51030, 51200,
584  51450, 51744, 51840, 52000, 52416, 52488, 52500, 52650, 52800, 52822,
585  52920, 53248, 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54880,
586  55000, 55296, 55440, 55566, 56000, 56160, 56250, 56320, 56448, 56700,
587  56862, 57024, 57330, 57344, 57600, 57624, 57750, 58212, 58240, 58320,
588  58500, 58800, 58968, 59136, 59400, 59904, 60000, 60368, 60480, 60750,
589  61152, 61236, 61250, 61440, 61600, 61740, 62208, 62370, 62400, 62426,
590  62500, 62720, 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512,
591  64680, 64800, 65000, 65520, 65536, 65610, 65856, 66000, 66150, 66528,
592  66560, 67200, 67228, 67392, 67500, 67584, 67914, 68040, 68250, 68600,
593  68750, 68796, 68992, 69120, 69300, 69888, 69984, 70000, 70200, 70400,
594  70560, 71280, 71344, 71442, 71680, 72000, 72030, 72576, 72800, 72900,
595  73500, 73710, 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264,
596  75460, 75600, 75816, 76032, 76440, 76800, 76832, 77000, 77616, 77760,
597  78000, 78400, 78624, 78732, 78750, 78848, 79200, 79380, 79872, 80000,
598  80190, 80262, 80640, 80850, 81000, 81250, 81536, 81648, 81900, 81920,
599  82320, 82500, 82944, 83160, 83200, 84000, 84240, 84480, 84672, 85050,
600  85536, 85750, 86016, 86240, 86400, 86436, 87318, 87360, 87480, 87500,
601  87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856,
602  90000, 90112, 90552, 90720, 91000, 91728, 91854, 92160, 92400, 92610,
603  93184, 93312, 93600, 93750, 94080, 94500, 94770, 95040, 95256, 95550,
604  96000, 96040, 96228, 96250, 96768, 97020, 97200, 97500, 98000, 98280,
605  98304, 98560, 98784, 99000, 99792, 99840 };
606 
607  if(s <= 0 || s > goodSizes[size-1])
608  return s;
609  return *std::upper_bound(goodSizes, goodSizes+size, s, std::less_equal<int>());
610 }
611 
612 template <int M>
613 struct FFTEmbedKernel
614 {
615  template <unsigned int N, class Real, class C, class Shape>
616  static void
617  exec(MultiArrayView<N, Real, C> & out, Shape const & kernelShape,
618  Shape & srcPoint, Shape & destPoint, bool copyIt)
619  {
620  for(srcPoint[M]=0; srcPoint[M]<kernelShape[M]; ++srcPoint[M])
621  {
622  if(srcPoint[M] < (kernelShape[M] + 1) / 2)
623  {
624  destPoint[M] = srcPoint[M];
625  }
626  else
627  {
628  destPoint[M] = srcPoint[M] + out.shape(M) - kernelShape[M];
629  copyIt = true;
630  }
631  FFTEmbedKernel<M-1>::exec(out, kernelShape, srcPoint, destPoint, copyIt);
632  }
633  }
634 };
635 
636 template <>
637 struct FFTEmbedKernel<0>
638 {
639  template <unsigned int N, class Real, class C, class Shape>
640  static void
641  exec(MultiArrayView<N, Real, C> & out, Shape const & kernelShape,
642  Shape & srcPoint, Shape & destPoint, bool copyIt)
643  {
644  for(srcPoint[0]=0; srcPoint[0]<kernelShape[0]; ++srcPoint[0])
645  {
646  if(srcPoint[0] < (kernelShape[0] + 1) / 2)
647  {
648  destPoint[0] = srcPoint[0];
649  }
650  else
651  {
652  destPoint[0] = srcPoint[0] + out.shape(0) - kernelShape[0];
653  copyIt = true;
654  }
655  if(copyIt)
656  {
657  out[destPoint] = out[srcPoint];
658  out[srcPoint] = 0.0;
659  }
660  }
661  }
662 };
663 
664 template <unsigned int N, class Real, class C1, class C2>
665 void
666 fftEmbedKernel(MultiArrayView<N, Real, C1> kernel,
667  MultiArrayView<N, Real, C2> out,
668  Real norm = 1.0)
669 {
670  typedef typename MultiArrayShape<N>::type Shape;
671 
672  MultiArrayView<N, Real, C2> kout = out.subarray(Shape(), kernel.shape());
673 
674  out.init(0.0);
675  kout = kernel;
676  if (norm != 1.0)
677  kout *= norm;
678  moveDCToUpperLeft(kout);
679 
680  Shape srcPoint, destPoint;
681  FFTEmbedKernel<(int)N-1>::exec(out, kernel.shape(), srcPoint, destPoint, false);
682 }
683 
684 template <unsigned int N, class Real, class C1, class C2>
685 void
686 fftEmbedArray(MultiArrayView<N, Real, C1> in,
687  MultiArrayView<N, Real, C2> out)
688 {
689  typedef typename MultiArrayShape<N>::type Shape;
690 
691  Shape diff = out.shape() - in.shape(),
692  leftDiff = div(diff, MultiArrayIndex(2)),
693  rightDiff = diff - leftDiff,
694  right = in.shape() + leftDiff;
695 
696  out.subarray(leftDiff, right) = in;
697 
698  typedef typename MultiArrayView<N, Real, C2>::traverser Traverser;
699  typedef MultiArrayNavigator<Traverser, N> Navigator;
700  typedef typename Navigator::iterator Iterator;
701 
702  for(unsigned int d = 0; d < N; ++d)
703  {
704  Navigator nav(out.traverser_begin(), out.shape(), d);
705 
706  for( ; nav.hasMore(); nav++ )
707  {
708  Iterator i = nav.begin();
709  for(int k=1; k<=leftDiff[d]; ++k)
710  i[leftDiff[d] - k] = i[leftDiff[d] + k];
711  for(int k=0; k<rightDiff[d]; ++k)
712  i[right[d] + k] = i[right[d] - k - 2];
713  }
714  }
715 }
716 
717 } // namespace detail
718 
719 template <class T, int N>
720 TinyVector<T, N>
721 fftwBestPaddedShape(TinyVector<T, N> shape)
722 {
723  for(unsigned int k=0; k<N; ++k)
724  shape[k] = detail::fftwPaddingSize(shape[k]);
725  return shape;
726 }
727 
728 template <class T, int N>
729 TinyVector<T, N>
730 fftwBestPaddedShapeR2C(TinyVector<T, N> shape)
731 {
732  shape[0] = detail::fftwEvenPaddingSize(shape[0]);
733  for(unsigned int k=1; k<N; ++k)
734  shape[k] = detail::fftwPaddingSize(shape[k]);
735  return shape;
736 }
737 
738 /** \brief Find frequency domain shape for a R2C Fourier transform.
739 
740  When a real valued array is transformed to the frequency domain, about half of the
741  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
742  transform</a> that doesn't compute and store the redundant coefficients. This function
743  computes the appropriate frequency domain shape for a given shape in the spatial domain.
744  It simply replaces <tt>shape[0]</tt> with <tt>shape[0] / 2 + 1</tt>.
745 
746  <b>\#include</b> <vigra/multi_fft.hxx>
747 */
748 template <class T, int N>
749 TinyVector<T, N>
751 {
752  shape[0] = shape[0] / 2 + 1;
753  return shape;
754 }
755 
756 template <class T, int N>
757 TinyVector<T, N>
758 fftwCorrespondingShapeC2R(TinyVector<T, N> shape, bool oddDimension0 = false)
759 {
760  shape[0] = oddDimension0
761  ? (shape[0] - 1) * 2 + 1
762  : (shape[0] - 1) * 2;
763  return shape;
764 }
765 
766 /********************************************************/
767 /* */
768 /* FFTWPlan */
769 /* */
770 /********************************************************/
771 
772 /** C++ wrapper for FFTW plans.
773 
774  The class encapsulates the calls to <tt>fftw_plan_dft_2d</tt>, <tt>fftw_execute</tt>, and
775  <tt>fftw_destroy_plan</tt> (and their <tt>float</tt> and <tt>long double</tt> counterparts)
776  in an easy-to-use interface.
777 
778  Usually, you use this class only indirectly via \ref fourierTransform()
779  and \ref fourierTransformInverse(). You only need this class if you want to have more control
780  about FFTW's planning process (by providing non-default planning flags) and/or want to re-use
781  plans for several transformations.
782 
783  <b> Usage:</b>
784 
785  <b>\#include</b> <vigra/multi_fft.hxx><br>
786  Namespace: vigra
787 
788  \code
789  // compute complex Fourier transform of a real image
790  MultiArray<2, double> src(Shape2(w, h));
791  MultiArray<2, FFTWComplex<double> > fourier(Shape2(w, h));
792 
793  // create an optimized plan by measuring the speed of several algorithm variants
794  FFTWPlan<2, double> plan(src, fourier, FFTW_MEASURE);
795 
796  plan.execute(src, fourier);
797  \endcode
798 */
799 template <unsigned int N, class Real = double>
800 class FFTWPlan
801 {
802  typedef ArrayVector<int> Shape;
803  typedef typename FFTWReal2Complex<Real>::plan_type PlanType;
804  typedef typename FFTWComplex<Real>::complex_type Complex;
805 
806  PlanType plan;
807  Shape shape, instrides, outstrides;
808  int sign;
809 
810  public:
811  /** \brief Create an empty plan.
812 
813  The plan can be initialized later by one of the init() functions.
814  */
816  : plan(0)
817  {}
818 
819  /** \brief Create a plan for a complex-to-complex transform.
820 
821  \arg SIGN must be <tt>FFTW_FORWARD</tt> or <tt>FFTW_BACKWARD</tt> according to the
822  desired transformation direction.
823  \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner
824  flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess
825  optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>.
826  */
827  template <class C1, class C2>
829  MultiArrayView<N, FFTWComplex<Real>, C2> out,
830  int SIGN, unsigned int planner_flags = FFTW_ESTIMATE)
831  : plan(0)
832  {
833  init(in, out, SIGN, planner_flags);
834  }
835 
836  /** \brief Create a plan for a real-to-complex transform.
837 
838  This always refers to a forward transform. The shape of the output determines
839  if a standard transform (when <tt>out.shape() == in.shape()</tt>) or an
840  <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C
841  transform</a> (when <tt>out.shape() == fftwCorrespondingShapeR2C(in.shape())</tt>) will be executed.
842 
843  \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner
844  flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess
845  optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>.
846  */
847  template <class C1, class C2>
849  MultiArrayView<N, FFTWComplex<Real>, C2> out,
850  unsigned int planner_flags = FFTW_ESTIMATE)
851  : plan(0)
852  {
853  init(in, out, planner_flags);
854  }
855 
856  /** \brief Create a plan for a complex-to-real transform.
857 
858  This always refers to a inverse transform. The shape of the input determines
859  if a standard transform (when <tt>in.shape() == out.shape()</tt>) or a
860  <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">C2R
861  transform</a> (when <tt>in.shape() == fftwCorrespondingShapeR2C(out.shape())</tt>) will be executed.
862 
863  \arg planner_flags must be a combination of the <a href="http://www.fftw.org/doc/Planner-Flags.html">planner
864  flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess
865  optimal algorithm settings or read them from pre-loaded <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>.
866  */
867  template <class C1, class C2>
870  unsigned int planner_flags = FFTW_ESTIMATE)
871  : plan(0)
872  {
873  init(in, out, planner_flags);
874  }
875 
876  /** \brief Copy constructor.
877  */
878  FFTWPlan(FFTWPlan const & other)
879  : plan(other.plan),
880  sign(other.sign)
881  {
882  FFTWPlan & o = const_cast<FFTWPlan &>(other);
883  shape.swap(o.shape);
884  instrides.swap(o.instrides);
885  outstrides.swap(o.outstrides);
886  o.plan = 0; // act like std::auto_ptr
887  }
888 
889  /** \brief Copy assigment.
890  */
891  FFTWPlan & operator=(FFTWPlan const & other)
892  {
893  if(this != &other)
894  {
895  FFTWPlan & o = const_cast<FFTWPlan &>(other);
896  plan = o.plan;
897  shape.swap(o.shape);
898  instrides.swap(o.instrides);
899  outstrides.swap(o.outstrides);
900  sign = o.sign;
901  o.plan = 0; // act like std::auto_ptr
902  }
903  return *this;
904  }
905 
906  /** \brief Destructor.
907  */
909  {
910  detail::fftwPlanDestroy(plan);
911  }
912 
913  /** \brief Init a complex-to-complex transform.
914 
915  See the constructor with the same signature for details.
916  */
917  template <class C1, class C2>
919  MultiArrayView<N, FFTWComplex<Real>, C2> out,
920  int SIGN, unsigned int planner_flags = FFTW_ESTIMATE)
921  {
922  vigra_precondition(in.strideOrdering() == out.strideOrdering(),
923  "FFTWPlan.init(): input and output must have the same stride ordering.");
924 
925  initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(),
926  SIGN, planner_flags);
927  }
928 
929  /** \brief Init a real-to-complex transform.
930 
931  See the constructor with the same signature for details.
932  */
933  template <class C1, class C2>
935  MultiArrayView<N, FFTWComplex<Real>, C2> out,
936  unsigned int planner_flags = FFTW_ESTIMATE)
937  {
938  vigra_precondition(in.strideOrdering() == out.strideOrdering(),
939  "FFTWPlan.init(): input and output must have the same stride ordering.");
940 
941  initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(),
942  FFTW_FORWARD, planner_flags);
943  }
944 
945  /** \brief Init a complex-to-real transform.
946 
947  See the constructor with the same signature for details.
948  */
949  template <class C1, class C2>
952  unsigned int planner_flags = FFTW_ESTIMATE)
953  {
954  vigra_precondition(in.strideOrdering() == out.strideOrdering(),
955  "FFTWPlan.init(): input and output must have the same stride ordering.");
956 
957  initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(),
958  FFTW_BACKWARD, planner_flags);
959  }
960 
961  /** \brief Execute a complex-to-complex transform.
962 
963  The array shapes must be the same as in the corresponding init function
964  or constructor. However, execute() can be called several times on
965  the same plan, even with different arrays, as long as they have the appropriate
966  shapes.
967  */
968  template <class C1, class C2>
970  MultiArrayView<N, FFTWComplex<Real>, C2> out) const
971  {
972  executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending());
973  }
974 
975  /** \brief Execute a real-to-complex transform.
976 
977  The array shapes must be the same as in the corresponding init function
978  or constructor. However, execute() can be called several times on
979  the same plan, even with different arrays, as long as they have the appropriate
980  shapes.
981  */
982  template <class C1, class C2>
984  MultiArrayView<N, FFTWComplex<Real>, C2> out) const
985  {
986  executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending());
987  }
988 
989  /** \brief Execute a complex-to-real transform.
990 
991  The array shapes must be the same as in the corresponding init function
992  or constructor. However, execute() can be called several times on
993  the same plan, even with different arrays, as long as they have the appropriate
994  shapes.
995  */
996  template <class C1, class C2>
998  MultiArrayView<N, Real, C2> out) const
999  {
1000  executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending());
1001  }
1002 
1003  private:
1004 
1005  template <class MI, class MO>
1006  void initImpl(MI ins, MO outs, int SIGN, unsigned int planner_flags);
1007 
1008  template <class MI, class MO>
1009  void executeImpl(MI ins, MO outs) const;
1010 
1011  void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> in,
1013  {
1014  vigra_precondition(in.shape() == out.shape(),
1015  "FFTWPlan.init(): input and output must have the same shape.");
1016  }
1017 
1018  void checkShapes(MultiArrayView<N, Real, StridedArrayTag> ins,
1019  MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> outs) const
1020  {
1021  for(int k=0; k<(int)N-1; ++k)
1022  vigra_precondition(ins.shape(k) == outs.shape(k),
1023  "FFTWPlan.init(): input and output must have matching shapes.");
1024  vigra_precondition(ins.shape(N-1) / 2 + 1 == outs.shape(N-1),
1025  "FFTWPlan.init(): input and output must have matching shapes.");
1026  }
1027 
1028  void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> ins,
1029  MultiArrayView<N, Real, StridedArrayTag> outs) const
1030  {
1031  for(int k=0; k<(int)N-1; ++k)
1032  vigra_precondition(ins.shape(k) == outs.shape(k),
1033  "FFTWPlan.init(): input and output must have matching shapes.");
1034  vigra_precondition(outs.shape(N-1) / 2 + 1 == ins.shape(N-1),
1035  "FFTWPlan.init(): input and output must have matching shapes.");
1036  }
1037 };
1038 
1039 template <unsigned int N, class Real>
1040 template <class MI, class MO>
1041 void
1042 FFTWPlan<N, Real>::initImpl(MI ins, MO outs, int SIGN, unsigned int planner_flags)
1043 {
1044  checkShapes(ins, outs);
1045 
1046  typename MultiArrayShape<N>::type logicalShape(SIGN == FFTW_FORWARD
1047  ? ins.shape()
1048  : outs.shape());
1049 
1050  Shape newShape(logicalShape.begin(), logicalShape.end()),
1051  newIStrides(ins.stride().begin(), ins.stride().end()),
1052  newOStrides(outs.stride().begin(), outs.stride().end()),
1053  itotal(ins.shape().begin(), ins.shape().end()),
1054  ototal(outs.shape().begin(), outs.shape().end());
1055 
1056  for(unsigned int j=1; j<N; ++j)
1057  {
1058  itotal[j] = ins.stride(j-1) / ins.stride(j);
1059  ototal[j] = outs.stride(j-1) / outs.stride(j);
1060  }
1061 
1062  PlanType newPlan = detail::fftwPlanCreate(N, newShape.begin(),
1063  ins.data(), itotal.begin(), ins.stride(N-1),
1064  outs.data(), ototal.begin(), outs.stride(N-1),
1065  SIGN, planner_flags);
1066  detail::fftwPlanDestroy(plan);
1067  plan = newPlan;
1068  shape.swap(newShape);
1069  instrides.swap(newIStrides);
1070  outstrides.swap(newOStrides);
1071  sign = SIGN;
1072 }
1073 
1074 template <unsigned int N, class Real>
1075 template <class MI, class MO>
1076 void FFTWPlan<N, Real>::executeImpl(MI ins, MO outs) const
1077 {
1078  vigra_precondition(plan != 0, "FFTWPlan::execute(): plan is NULL.");
1079 
1080  typename MultiArrayShape<N>::type lshape(sign == FFTW_FORWARD
1081  ? ins.shape()
1082  : outs.shape());
1083 
1084  vigra_precondition((lshape == TinyVectorView<int, N>(shape.data())),
1085  "FFTWPlan::execute(): shape mismatch between plan and data.");
1086  vigra_precondition((ins.stride() == TinyVectorView<int, N>(instrides.data())),
1087  "FFTWPlan::execute(): strides mismatch between plan and input data.");
1088  vigra_precondition((outs.stride() == TinyVectorView<int, N>(outstrides.data())),
1089  "FFTWPlan::execute(): strides mismatch between plan and output data.");
1090 
1091  detail::fftwPlanExecute(plan, ins.data(), outs.data());
1092 
1093  typedef typename MO::value_type V;
1094  if(sign == FFTW_BACKWARD)
1095  outs *= V(1.0) / Real(outs.size());
1096 }
1097 
1098 /********************************************************/
1099 /* */
1100 /* FFTWConvolvePlan */
1101 /* */
1102 /********************************************************/
1103 
1104 /** C++ wrapper for a pair of FFTW plans used to perform FFT-based convolution.
1105 
1106  The class encapsulates the calls to <tt>fftw_plan_dft_2d</tt>, <tt>fftw_execute</tt>, and
1107  <tt>fftw_destroy_plan</tt> (and their <tt>float</tt> and <tt>long double</tt> counterparts)
1108  in an easy-to-use interface. It always creates a pair of plans, one for the forward and one
1109  for the inverse transform required for convolution.
1110 
1111  Usually, you use this class only indirectly via \ref convolveFFT() and its variants.
1112  You only need this class if you want to have more control about FFTW's planning process
1113  (by providing non-default planning flags) and/or want to re-use plans for several convolutions.
1114 
1115  <b> Usage:</b>
1116 
1117  <b>\#include</b> <vigra/multi_fft.hxx><br>
1118  Namespace: vigra
1119 
1120  \code
1121  // convolve a real array with a real kernel
1122  MultiArray<2, double> src(Shape2(w, h)), dest(Shape2(w, h));
1123 
1124  MultiArray<2, double> spatial_kernel(Shape2(9, 9));
1125  Gaussian<double> gauss(1.0);
1126 
1127  for(int y=0; y<9; ++y)
1128  for(int x=0; x<9; ++x)
1129  spatial_kernel(x, y) = gauss(x-4.0)*gauss(y-4.0);
1130 
1131  // create an optimized plan by measuring the speed of several algorithm variants
1132  FFTWConvolvePlan<2, double> plan(src, spatial_kernel, dest, FFTW_MEASURE);
1133 
1134  plan.execute(src, spatial_kernel, dest);
1135  \endcode
1136 */
1137 template <unsigned int N, class Real>
1139 {
1140  typedef FFTWComplex<Real> Complex;
1143 
1144  FFTWPlan<N, Real> forward_plan, backward_plan;
1145  RArray realArray, realKernel;
1146  CArray fourierArray, fourierKernel;
1147  bool useFourierKernel;
1148 
1149  public:
1150 
1151  typedef typename MultiArrayShape<N>::type Shape;
1152 
1153  /** \brief Create an empty plan.
1154 
1155  The plan can be initialized later by one of the init() functions.
1156  */
1158  : useFourierKernel(false)
1159  {}
1160 
1161  /** \brief Create a plan to convolve a real array with a real kernel.
1162 
1163  The kernel must be defined in the spatial domain.
1164  See \ref convolveFFT() for detailed information on required shapes and internal padding.
1165 
1166  \arg planner_flags must be a combination of the
1167  <a href="http://www.fftw.org/doc/Planner-Flags.html">planner
1168  flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess
1169  optimal algorithm settings or read them from pre-loaded
1170  <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>.
1171  */
1172  template <class C1, class C2, class C3>
1176  unsigned int planner_flags = FFTW_ESTIMATE)
1177  : useFourierKernel(false)
1178  {
1179  init(in, kernel, out, planner_flags);
1180  }
1181 
1182  /** \brief Create a plan to convolve a real array with a complex kernel.
1183 
1184  The kernel must be defined in the Fourier domain, using the half-space format.
1185  See \ref convolveFFT() for detailed information on required shapes and internal padding.
1186 
1187  \arg planner_flags must be a combination of the
1188  <a href="http://www.fftw.org/doc/Planner-Flags.html">planner
1189  flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess
1190  optimal algorithm settings or read them from pre-loaded
1191  <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>.
1192  */
1193  template <class C1, class C2, class C3>
1195  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1197  unsigned int planner_flags = FFTW_ESTIMATE)
1198  : useFourierKernel(true)
1199  {
1200  init(in, kernel, out, planner_flags);
1201  }
1202 
1203  /** \brief Create a plan to convolve a complex array with a complex kernel.
1204 
1205  See \ref convolveFFT() for detailed information on required shapes and internal padding.
1206 
1207  \arg fourierDomainKernel determines if the kernel is defined in the spatial or
1208  Fourier domain.
1209  \arg planner_flags must be a combination of the
1210  <a href="http://www.fftw.org/doc/Planner-Flags.html">planner
1211  flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess
1212  optimal algorithm settings or read them from pre-loaded
1213  <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>.
1214  */
1215  template <class C1, class C2, class C3>
1217  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1218  MultiArrayView<N, FFTWComplex<Real>, C3> out,
1219  bool fourierDomainKernel,
1220  unsigned int planner_flags = FFTW_ESTIMATE)
1221  {
1222  init(in, kernel, out, fourierDomainKernel, planner_flags);
1223  }
1224 
1225 
1226  /** \brief Create a plan from just the shape information.
1227 
1228  See \ref convolveFFT() for detailed information on required shapes and internal padding.
1229 
1230  \arg fourierDomainKernel determines if the kernel is defined in the spatial or
1231  Fourier domain.
1232  \arg planner_flags must be a combination of the
1233  <a href="http://www.fftw.org/doc/Planner-Flags.html">planner
1234  flags</a> defined by the FFTW library. The default <tt>FFTW_ESTIMATE</tt> will guess
1235  optimal algorithm settings or read them from pre-loaded
1236  <a href="http://www.fftw.org/doc/Wisdom.html">"wisdom"</a>.
1237  */
1238  template <class C1, class C2, class C3>
1239  FFTWConvolvePlan(Shape inOut, Shape kernel,
1240  bool useFourierKernel = false,
1241  unsigned int planner_flags = FFTW_ESTIMATE)
1242  {
1243  if(useFourierKernel)
1244  init(inOut, kernel, planner_flags);
1245  else
1246  initFourierKernel(inOut, kernel, planner_flags);
1247  }
1248 
1249  /** \brief Init a plan to convolve a real array with a real kernel.
1250 
1251  See the constructor with the same signature for details.
1252  */
1253  template <class C1, class C2, class C3>
1257  unsigned int planner_flags = FFTW_ESTIMATE)
1258  {
1259  vigra_precondition(in.shape() == out.shape(),
1260  "FFTWConvolvePlan::init(): input and output must have the same shape.");
1261  init(in.shape(), kernel.shape(), planner_flags);
1262  }
1263 
1264  /** \brief Init a plan to convolve a real array with a complex kernel.
1265 
1266  See the constructor with the same signature for details.
1267  */
1268  template <class C1, class C2, class C3>
1270  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1272  unsigned int planner_flags = FFTW_ESTIMATE)
1273  {
1274  vigra_precondition(in.shape() == out.shape(),
1275  "FFTWConvolvePlan::init(): input and output must have the same shape.");
1276  initFourierKernel(in.shape(), kernel.shape(), planner_flags);
1277  }
1278 
1279  /** \brief Init a plan to convolve a complex array with a complex kernel.
1280 
1281  See the constructor with the same signature for details.
1282  */
1283  template <class C1, class C2, class C3>
1285  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1286  MultiArrayView<N, FFTWComplex<Real>, C3> out,
1287  bool fourierDomainKernel,
1288  unsigned int planner_flags = FFTW_ESTIMATE)
1289  {
1290  vigra_precondition(in.shape() == out.shape(),
1291  "FFTWConvolvePlan::init(): input and output must have the same shape.");
1292  useFourierKernel = fourierDomainKernel;
1293  initComplex(in.shape(), kernel.shape(), planner_flags);
1294  }
1295 
1296  /** \brief Init a plan to convolve a real array with a sequence of kernels.
1297 
1298  The kernels can be either real or complex. The sequences \a kernels and \a outs
1299  must have the same length. See the corresponding constructors
1300  for single kernels for details.
1301  */
1302  template <class C1, class KernelIterator, class OutIterator>
1304  KernelIterator kernels, KernelIterator kernelsEnd,
1305  OutIterator outs, unsigned int planner_flags = FFTW_ESTIMATE)
1306  {
1307  typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1308  typedef typename KernelArray::value_type KernelValue;
1309  typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1310  typedef typename OutArray::value_type OutValue;
1311 
1312  bool realKernel = IsSameType<KernelValue, Real>::value;
1313  bool fourierKernel = IsSameType<KernelValue, Complex>::value;
1314 
1315  vigra_precondition(realKernel || fourierKernel,
1316  "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
1317  vigra_precondition((IsSameType<OutValue, Real>::value),
1318  "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
1319 
1320  if(realKernel)
1321  {
1322  initMany(in.shape(), checkShapes(in.shape(), kernels, kernelsEnd, outs),
1323  planner_flags);
1324  }
1325  else
1326  {
1327  initFourierKernelMany(in.shape(),
1328  checkShapesFourier(in.shape(), kernels, kernelsEnd, outs),
1329  planner_flags);
1330  }
1331  }
1332 
1333  /** \brief Init a plan to convolve a complex array with a sequence of kernels.
1334 
1335  The kernels must be complex as well. The sequences \a kernels and \a outs
1336  must have the same length. See the corresponding constructors
1337  for single kernels for details.
1338  */
1339  template <class C1, class KernelIterator, class OutIterator>
1341  KernelIterator kernels, KernelIterator kernelsEnd,
1342  OutIterator outs,
1343  bool fourierDomainKernels,
1344  unsigned int planner_flags = FFTW_ESTIMATE)
1345  {
1346  typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1347  typedef typename KernelArray::value_type KernelValue;
1348  typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1349  typedef typename OutArray::value_type OutValue;
1350 
1351  vigra_precondition((IsSameType<KernelValue, Complex>::value),
1352  "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
1353  vigra_precondition((IsSameType<OutValue, Complex>::value),
1354  "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
1355 
1356  useFourierKernel = fourierDomainKernels;
1357 
1358  Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs);
1359 
1360  CArray newFourierArray(paddedShape), newFourierKernel(paddedShape);
1361 
1362  FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags);
1363  FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags);
1364 
1365  forward_plan = fplan;
1366  backward_plan = bplan;
1367  fourierArray.swap(newFourierArray);
1368  fourierKernel.swap(newFourierKernel);
1369  }
1370 
1371  void init(Shape inOut, Shape kernel,
1372  unsigned int planner_flags = FFTW_ESTIMATE);
1373 
1374  void initFourierKernel(Shape inOut, Shape kernel,
1375  unsigned int planner_flags = FFTW_ESTIMATE);
1376 
1377  void initComplex(Shape inOut, Shape kernel,
1378  unsigned int planner_flags = FFTW_ESTIMATE);
1379 
1380  void initMany(Shape inOut, Shape maxKernel,
1381  unsigned int planner_flags = FFTW_ESTIMATE)
1382  {
1383  init(inOut, maxKernel, planner_flags);
1384  }
1385 
1386  void initFourierKernelMany(Shape inOut, Shape kernels,
1387  unsigned int planner_flags = FFTW_ESTIMATE)
1388  {
1389  initFourierKernel(inOut, kernels, planner_flags);
1390  }
1391 
1392  /** \brief Execute a plan to convolve a real array with a real kernel.
1393 
1394  The array shapes must be the same as in the corresponding init function
1395  or constructor. However, execute() can be called several times on
1396  the same plan, even with different arrays, as long as they have the appropriate
1397  shapes.
1398  */
1399  template <class C1, class C2, class C3>
1400  void execute(MultiArrayView<N, Real, C1> in,
1401  MultiArrayView<N, Real, C2> kernel,
1402  MultiArrayView<N, Real, C3> out);
1403 
1404  /** \brief Execute a plan to convolve a real array with a complex kernel.
1405 
1406  The array shapes must be the same as in the corresponding init function
1407  or constructor. However, execute() can be called several times on
1408  the same plan, even with different arrays, as long as they have the appropriate
1409  shapes.
1410  */
1411  template <class C1, class C2, class C3>
1412  void execute(MultiArrayView<N, Real, C1> in,
1413  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1414  MultiArrayView<N, Real, C3> out);
1415 
1416  /** \brief Execute a plan to convolve a complex array with a complex kernel.
1417 
1418  The array shapes must be the same as in the corresponding init function
1419  or constructor. However, execute() can be called several times on
1420  the same plan, even with different arrays, as long as they have the appropriate
1421  shapes.
1422  */
1423  template <class C1, class C2, class C3>
1424  void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1425  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1426  MultiArrayView<N, FFTWComplex<Real>, C3> out);
1427 
1428 
1429  /** \brief Execute a plan to convolve a complex array with a sequence of kernels.
1430 
1431  The array shapes must be the same as in the corresponding init function
1432  or constructor. However, executeMany() can be called several times on
1433  the same plan, even with different arrays, as long as they have the appropriate
1434  shapes.
1435  */
1436  template <class C1, class KernelIterator, class OutIterator>
1437  void executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1438  KernelIterator kernels, KernelIterator kernelsEnd,
1439  OutIterator outs);
1440 
1441  /** \brief Execute a plan to convolve a real array with a sequence of kernels.
1442 
1443  The array shapes must be the same as in the corresponding init function
1444  or constructor. However, executeMany() can be called several times on
1445  the same plan, even with different arrays, as long as they have the appropriate
1446  shapes.
1447  */
1448  template <class C1, class KernelIterator, class OutIterator>
1450  KernelIterator kernels, KernelIterator kernelsEnd,
1451  OutIterator outs)
1452  {
1453  typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1454  typedef typename KernelArray::value_type KernelValue;
1455  typedef typename IsSameType<KernelValue, Complex>::type UseFourierKernel;
1456  typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1457  typedef typename OutArray::value_type OutValue;
1458 
1459  bool realKernel = IsSameType<KernelValue, Real>::value;
1460  bool fourierKernel = IsSameType<KernelValue, Complex>::value;
1461 
1462  vigra_precondition(realKernel || fourierKernel,
1463  "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type.");
1464  vigra_precondition((IsSameType<OutValue, Real>::value),
1465  "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type.");
1466 
1467  executeManyImpl(in, kernels, kernelsEnd, outs, UseFourierKernel());
1468  }
1469 
1470  private:
1471 
1472  template <class KernelIterator, class OutIterator>
1473  Shape checkShapes(Shape in,
1474  KernelIterator kernels, KernelIterator kernelsEnd,
1475  OutIterator outs);
1476 
1477  template <class KernelIterator, class OutIterator>
1478  Shape checkShapesFourier(Shape in,
1479  KernelIterator kernels, KernelIterator kernelsEnd,
1480  OutIterator outs);
1481 
1482  template <class KernelIterator, class OutIterator>
1483  Shape checkShapesComplex(Shape in,
1484  KernelIterator kernels, KernelIterator kernelsEnd,
1485  OutIterator outs);
1486 
1487  template <class C1, class KernelIterator, class OutIterator>
1488  void
1489  executeManyImpl(MultiArrayView<N, Real, C1> in,
1490  KernelIterator kernels, KernelIterator kernelsEnd,
1491  OutIterator outs, VigraFalseType /* useFourierKernel*/);
1492 
1493  template <class C1, class KernelIterator, class OutIterator>
1494  void
1495  executeManyImpl(MultiArrayView<N, Real, C1> in,
1496  KernelIterator kernels, KernelIterator kernelsEnd,
1497  OutIterator outs, VigraTrueType /* useFourierKernel*/);
1498 
1499 };
1500 
1501 template <unsigned int N, class Real>
1502 void
1503 FFTWConvolvePlan<N, Real>::init(Shape in, Shape kernel,
1504  unsigned int planner_flags)
1505 {
1506  Shape paddedShape = fftwBestPaddedShapeR2C(in + kernel - Shape(1)),
1507  complexShape = fftwCorrespondingShapeR2C(paddedShape);
1508 
1509  CArray newFourierArray(complexShape), newFourierKernel(complexShape);
1510 
1511  Shape realStrides = 2*newFourierArray.stride();
1512  realStrides[0] = 1;
1513  RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data());
1514  RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data());
1515 
1516  FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags);
1517  FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags);
1518 
1519  forward_plan = fplan;
1520  backward_plan = bplan;
1521  realArray = newRealArray;
1522  realKernel = newRealKernel;
1523  fourierArray.swap(newFourierArray);
1524  fourierKernel.swap(newFourierKernel);
1525  useFourierKernel = false;
1526 }
1527 
1528 template <unsigned int N, class Real>
1529 void
1530 FFTWConvolvePlan<N, Real>::initFourierKernel(Shape in, Shape kernel,
1531  unsigned int planner_flags)
1532 {
1533  Shape complexShape = kernel,
1534  paddedShape = fftwCorrespondingShapeC2R(complexShape);
1535 
1536  for(unsigned int k=0; k<N; ++k)
1537  vigra_precondition(in[k] <= paddedShape[k],
1538  "FFTWConvolvePlan::init(): kernel too small for given input.");
1539 
1540  CArray newFourierArray(complexShape), newFourierKernel(complexShape);
1541 
1542  Shape realStrides = 2*newFourierArray.stride();
1543  realStrides[0] = 1;
1544  RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data());
1545  RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data());
1546 
1547  FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags);
1548  FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags);
1549 
1550  forward_plan = fplan;
1551  backward_plan = bplan;
1552  realArray = newRealArray;
1553  realKernel = newRealKernel;
1554  fourierArray.swap(newFourierArray);
1555  fourierKernel.swap(newFourierKernel);
1556  useFourierKernel = true;
1557 }
1558 
1559 template <unsigned int N, class Real>
1560 void
1561 FFTWConvolvePlan<N, Real>::initComplex(Shape in, Shape kernel,
1562  unsigned int planner_flags)
1563 {
1564  Shape paddedShape;
1565 
1566  if(useFourierKernel)
1567  {
1568  for(unsigned int k=0; k<N; ++k)
1569  vigra_precondition(in[k] <= kernel[k],
1570  "FFTWConvolvePlan::init(): kernel too small for given input.");
1571 
1572  paddedShape = kernel;
1573  }
1574  else
1575  {
1576  paddedShape = fftwBestPaddedShape(in + kernel - Shape(1));
1577  }
1578 
1579  CArray newFourierArray(paddedShape), newFourierKernel(paddedShape);
1580 
1581  FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags);
1582  FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags);
1583 
1584  forward_plan = fplan;
1585  backward_plan = bplan;
1586  fourierArray.swap(newFourierArray);
1587  fourierKernel.swap(newFourierKernel);
1588 }
1589 
1590 #ifndef DOXYGEN // doxygen documents these functions as free functions
1591 
1592 template <unsigned int N, class Real>
1593 template <class C1, class C2, class C3>
1594 void
1595 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, Real, C1> in,
1596  MultiArrayView<N, Real, C2> kernel,
1597  MultiArrayView<N, Real, C3> out)
1598 {
1599  vigra_precondition(!useFourierKernel,
1600  "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel.");
1601 
1602  vigra_precondition(in.shape() == out.shape(),
1603  "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1604 
1605  Shape paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernel.shape() - Shape(1)),
1606  diff = paddedShape - in.shape(),
1607  left = div(diff, MultiArrayIndex(2)),
1608  right = in.shape() + left;
1609 
1610  vigra_precondition(paddedShape == realArray.shape(),
1611  "FFTWConvolvePlan::execute(): shape mismatch between input and plan.");
1612 
1613  detail::fftEmbedArray(in, realArray);
1614  forward_plan.execute(realArray, fourierArray);
1615 
1616  detail::fftEmbedKernel(kernel, realKernel);
1617  forward_plan.execute(realKernel, fourierKernel);
1618 
1619  fourierArray *= fourierKernel;
1620 
1621  backward_plan.execute(fourierArray, realArray);
1622 
1623  out = realArray.subarray(left, right);
1624 }
1625 
1626 template <unsigned int N, class Real>
1627 template <class C1, class C2, class C3>
1628 void
1629 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, Real, C1> in,
1630  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1631  MultiArrayView<N, Real, C3> out)
1632 {
1633  vigra_precondition(useFourierKernel,
1634  "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel.");
1635 
1636  vigra_precondition(in.shape() == out.shape(),
1637  "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1638 
1639  vigra_precondition(kernel.shape() == fourierArray.shape(),
1640  "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan.");
1641 
1642  Shape paddedShape = fftwCorrespondingShapeC2R(kernel.shape(), odd(in.shape(0))),
1643  diff = paddedShape - in.shape(),
1644  left = div(diff, MultiArrayIndex(2)),
1645  right = in.shape() + left;
1646 
1647  vigra_precondition(paddedShape == realArray.shape(),
1648  "FFTWConvolvePlan::execute(): shape mismatch between input and plan.");
1649 
1650  detail::fftEmbedArray(in, realArray);
1651  forward_plan.execute(realArray, fourierArray);
1652 
1653  fourierKernel = kernel;
1654  moveDCToHalfspaceUpperLeft(fourierKernel);
1655 
1656  fourierArray *= fourierKernel;
1657 
1658  backward_plan.execute(fourierArray, realArray);
1659 
1660  out = realArray.subarray(left, right);
1661 }
1662 
1663 template <unsigned int N, class Real>
1664 template <class C1, class C2, class C3>
1665 void
1666 FFTWConvolvePlan<N, Real>::execute(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1667  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1668  MultiArrayView<N, FFTWComplex<Real>, C3> out)
1669 {
1670  vigra_precondition(in.shape() == out.shape(),
1671  "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1672 
1673  Shape paddedShape = fourierArray.shape(),
1674  diff = paddedShape - in.shape(),
1675  left = div(diff, MultiArrayIndex(2)),
1676  right = in.shape() + left;
1677 
1678  if(useFourierKernel)
1679  {
1680  vigra_precondition(kernel.shape() == fourierArray.shape(),
1681  "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan.");
1682 
1683  fourierKernel = kernel;
1684  moveDCToUpperLeft(fourierKernel);
1685  }
1686  else
1687  {
1688  detail::fftEmbedKernel(kernel, fourierKernel);
1689  forward_plan.execute(fourierKernel, fourierKernel);
1690  }
1691 
1692  detail::fftEmbedArray(in, fourierArray);
1693  forward_plan.execute(fourierArray, fourierArray);
1694 
1695  fourierArray *= fourierKernel;
1696 
1697  backward_plan.execute(fourierArray, fourierArray);
1698 
1699  out = fourierArray.subarray(left, right);
1700 }
1701 
1702 template <unsigned int N, class Real>
1703 template <class C1, class KernelIterator, class OutIterator>
1704 void
1705 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
1706  KernelIterator kernels, KernelIterator kernelsEnd,
1707  OutIterator outs, VigraFalseType /*useFourierKernel*/)
1708 {
1709  vigra_precondition(!useFourierKernel,
1710  "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel.");
1711 
1712  Shape kernelMax = checkShapes(in.shape(), kernels, kernelsEnd, outs),
1713  paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernelMax - Shape(1)),
1714  diff = paddedShape - in.shape(),
1715  left = div(diff, MultiArrayIndex(2)),
1716  right = in.shape() + left;
1717 
1718  vigra_precondition(paddedShape == realArray.shape(),
1719  "FFTWConvolvePlan::executeMany(): shape mismatch between input and plan.");
1720 
1721  detail::fftEmbedArray(in, realArray);
1722  forward_plan.execute(realArray, fourierArray);
1723 
1724  for(; kernels != kernelsEnd; ++kernels, ++outs)
1725  {
1726  detail::fftEmbedKernel(*kernels, realKernel);
1727  forward_plan.execute(realKernel, fourierKernel);
1728 
1729  fourierKernel *= fourierArray;
1730 
1731  backward_plan.execute(fourierKernel, realKernel);
1732 
1733  *outs = realKernel.subarray(left, right);
1734  }
1735 }
1736 
1737 template <unsigned int N, class Real>
1738 template <class C1, class KernelIterator, class OutIterator>
1739 void
1740 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
1741  KernelIterator kernels, KernelIterator kernelsEnd,
1742  OutIterator outs, VigraTrueType /*useFourierKernel*/)
1743 {
1744  vigra_precondition(useFourierKernel,
1745  "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel.");
1746 
1747  Shape complexShape = checkShapesFourier(in.shape(), kernels, kernelsEnd, outs),
1748  paddedShape = fftwCorrespondingShapeC2R(complexShape, odd(in.shape(0))),
1749  diff = paddedShape - in.shape(),
1750  left = div(diff, MultiArrayIndex(2)),
1751  right = in.shape() + left;
1752 
1753  vigra_precondition(complexShape == fourierArray.shape(),
1754  "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between kernels and plan.");
1755 
1756  vigra_precondition(paddedShape == realArray.shape(),
1757  "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between input and plan.");
1758 
1759  detail::fftEmbedArray(in, realArray);
1760  forward_plan.execute(realArray, fourierArray);
1761 
1762  for(; kernels != kernelsEnd; ++kernels, ++outs)
1763  {
1764  fourierKernel = *kernels;
1765  moveDCToHalfspaceUpperLeft(fourierKernel);
1766  fourierKernel *= fourierArray;
1767 
1768  backward_plan.execute(fourierKernel, realKernel);
1769 
1770  *outs = realKernel.subarray(left, right);
1771  }
1772 }
1773 
1774 template <unsigned int N, class Real>
1775 template <class C1, class KernelIterator, class OutIterator>
1776 void
1777 FFTWConvolvePlan<N, Real>::executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1778  KernelIterator kernels, KernelIterator kernelsEnd,
1779  OutIterator outs)
1780 {
1781  typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1782  typedef typename KernelArray::value_type KernelValue;
1783  typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1784  typedef typename OutArray::value_type OutValue;
1785 
1786  vigra_precondition((IsSameType<KernelValue, Complex>::value),
1787  "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type.");
1788  vigra_precondition((IsSameType<OutValue, Complex>::value),
1789  "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type.");
1790 
1791  Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs),
1792  diff = paddedShape - in.shape(),
1793  left = div(diff, MultiArrayIndex(2)),
1794  right = in.shape() + left;
1795 
1796  detail::fftEmbedArray(in, fourierArray);
1797  forward_plan.execute(fourierArray, fourierArray);
1798 
1799  for(; kernels != kernelsEnd; ++kernels, ++outs)
1800  {
1801  if(useFourierKernel)
1802  {
1803  fourierKernel = *kernels;
1804  moveDCToUpperLeft(fourierKernel);
1805  }
1806  else
1807  {
1808  detail::fftEmbedKernel(*kernels, fourierKernel);
1809  forward_plan.execute(fourierKernel, fourierKernel);
1810  }
1811 
1812  fourierKernel *= fourierArray;
1813 
1814  backward_plan.execute(fourierKernel, fourierKernel);
1815 
1816  *outs = fourierKernel.subarray(left, right);
1817  }
1818 }
1819 
1820 #endif // DOXYGEN
1821 
1822 template <unsigned int N, class Real>
1823 template <class KernelIterator, class OutIterator>
1824 typename FFTWConvolvePlan<N, Real>::Shape
1825 FFTWConvolvePlan<N, Real>::checkShapes(Shape in,
1826  KernelIterator kernels, KernelIterator kernelsEnd,
1827  OutIterator outs)
1828 {
1829  vigra_precondition(kernels != kernelsEnd,
1830  "FFTWConvolvePlan::checkShapes(): empty kernel sequence.");
1831 
1832  Shape kernelMax;
1833  for(; kernels != kernelsEnd; ++kernels, ++outs)
1834  {
1835  vigra_precondition(in == outs->shape(),
1836  "FFTWConvolvePlan::checkShapes(): shape mismatch between input and (one) output.");
1837  kernelMax = max(kernelMax, kernels->shape());
1838  }
1839  vigra_precondition(prod(kernelMax) > 0,
1840  "FFTWConvolvePlan::checkShapes(): all kernels have size 0.");
1841  return kernelMax;
1842 }
1843 
1844 template <unsigned int N, class Real>
1845 template <class KernelIterator, class OutIterator>
1846 typename FFTWConvolvePlan<N, Real>::Shape
1847 FFTWConvolvePlan<N, Real>::checkShapesFourier(Shape in,
1848  KernelIterator kernels, KernelIterator kernelsEnd,
1849  OutIterator outs)
1850 {
1851  vigra_precondition(kernels != kernelsEnd,
1852  "FFTWConvolvePlan::checkShapesFourier(): empty kernel sequence.");
1853 
1854  Shape complexShape = kernels->shape(),
1855  paddedShape = fftwCorrespondingShapeC2R(complexShape);
1856 
1857  for(unsigned int k=0; k<N; ++k)
1858  vigra_precondition(in[k] <= paddedShape[k],
1859  "FFTWConvolvePlan::checkShapesFourier(): kernels too small for given input.");
1860 
1861  for(; kernels != kernelsEnd; ++kernels, ++outs)
1862  {
1863  vigra_precondition(in == outs->shape(),
1864  "FFTWConvolvePlan::checkShapesFourier(): shape mismatch between input and (one) output.");
1865  vigra_precondition(complexShape == kernels->shape(),
1866  "FFTWConvolvePlan::checkShapesFourier(): all kernels must have the same size.");
1867  }
1868  return complexShape;
1869 }
1870 
1871 template <unsigned int N, class Real>
1872 template <class KernelIterator, class OutIterator>
1873 typename FFTWConvolvePlan<N, Real>::Shape
1874 FFTWConvolvePlan<N, Real>::checkShapesComplex(Shape in,
1875  KernelIterator kernels, KernelIterator kernelsEnd,
1876  OutIterator outs)
1877 {
1878  vigra_precondition(kernels != kernelsEnd,
1879  "FFTWConvolvePlan::checkShapesComplex(): empty kernel sequence.");
1880 
1881  Shape kernelShape = kernels->shape();
1882  for(; kernels != kernelsEnd; ++kernels, ++outs)
1883  {
1884  vigra_precondition(in == outs->shape(),
1885  "FFTWConvolvePlan::checkShapesComplex(): shape mismatch between input and (one) output.");
1886  if(useFourierKernel)
1887  {
1888  vigra_precondition(kernelShape == kernels->shape(),
1889  "FFTWConvolvePlan::checkShapesComplex(): Fourier domain kernels must have identical size.");
1890  }
1891  else
1892  {
1893  kernelShape = max(kernelShape, kernels->shape());
1894  }
1895  }
1896  vigra_precondition(prod(kernelShape) > 0,
1897  "FFTWConvolvePlan::checkShapesComplex(): all kernels have size 0.");
1898 
1899  if(useFourierKernel)
1900  {
1901  for(unsigned int k=0; k<N; ++k)
1902  vigra_precondition(in[k] <= kernelShape[k],
1903  "FFTWConvolvePlan::checkShapesComplex(): kernels too small for given input.");
1904  return kernelShape;
1905  }
1906  else
1907  {
1908  return fftwBestPaddedShape(in + kernelShape - Shape(1));
1909  }
1910 }
1911 
1912 
1913 /********************************************************/
1914 /* */
1915 /* fourierTransform */
1916 /* */
1917 /********************************************************/
1918 
1919 template <unsigned int N, class Real, class C1, class C2>
1920 inline void
1921 fourierTransform(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1922  MultiArrayView<N, FFTWComplex<Real>, C2> out)
1923 {
1924  FFTWPlan<N, Real>(in, out, FFTW_FORWARD).execute(in, out);
1925 }
1926 
1927 template <unsigned int N, class Real, class C1, class C2>
1928 inline void
1929 fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1930  MultiArrayView<N, FFTWComplex<Real>, C2> out)
1931 {
1932  FFTWPlan<N, Real>(in, out, FFTW_BACKWARD).execute(in, out);
1933 }
1934 
1935 template <unsigned int N, class Real, class C1, class C2>
1936 void
1937 fourierTransform(MultiArrayView<N, Real, C1> in,
1938  MultiArrayView<N, FFTWComplex<Real>, C2> out)
1939 {
1940  if(in.shape() == out.shape())
1941  {
1942  // copy the input array into the output and then perform an in-place FFT
1943  out = in;
1944  FFTWPlan<N, Real>(out, out, FFTW_FORWARD).execute(out, out);
1945  }
1946  else if(out.shape() == fftwCorrespondingShapeR2C(in.shape()))
1947  {
1948  FFTWPlan<N, Real>(in, out).execute(in, out);
1949  }
1950  else
1951  vigra_precondition(false,
1952  "fourierTransform(): shape mismatch between input and output.");
1953 }
1954 
1955 template <unsigned int N, class Real, class C1, class C2>
1956 void
1957 fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1958  MultiArrayView<N, Real, C2> out)
1959 {
1960  vigra_precondition(in.shape() == fftwCorrespondingShapeR2C(out.shape()),
1961  "fourierTransformInverse(): shape mismatch between input and output.");
1962  FFTWPlan<N, Real>(in, out).execute(in, out);
1963 }
1964 
1965 //@}
1966 
1967 /** \addtogroup MultiArrayConvolutionFilters
1968 */
1969 //@{
1970 
1971 /********************************************************/
1972 /* */
1973 /* convolveFFT */
1974 /* */
1975 /********************************************************/
1976 
1977 /** \brief Convolve an array with a kernel by means of the Fourier transform.
1978 
1979  Thanks to the convolution theorem of Fourier theory, a convolution in the spatial domain
1980  is equivalent to a multiplication in the frequency domain. Thus, for certain kernels
1981  (especially large, non-separable ones), it is advantageous to perform the convolution by first
1982  transforming both array and kernel to the frequency domain, multiplying the frequency
1983  representations, and transforming the result back into the spatial domain.
1984  Some kernels have a much simpler definition in the frequency domain, so that they are readily
1985  computed there directly, avoiding Fourier transformation of those kernels.
1986 
1987  The following functions implement various variants of FFT-based convolution:
1988 
1989  <DL>
1990  <DT><b>convolveFFT</b><DD> Convolve a real-valued input array with a kernel such that the
1991  result is also real-valued. That is, the kernel is either provided
1992  as a real-valued array in the spatial domain, or as a
1993  complex-valued array in the Fourier domain, using the half-space format
1994  of the R2C Fourier transform (see below).
1995  <DT><b>convolveFFTMany</b><DD> Like <tt>convolveFFT</tt>, but you may provide many kernels at once
1996  (using an iterator pair specifying the kernel sequence).
1997  This has the advantage that the forward transform of the input array needs
1998  to be executed only once.
1999  <DT><b>convolveFFTComplex</b><DD> Convolve a complex-valued input array with a complex-valued kernel,
2000  resulting in a complex-valued output array. An additional flag is used to
2001  specify whether the kernel is defined in the spatial or frequency domain.
2002  <DT><b>convolveFFTComplexMany</b><DD> Like <tt>convolveFFTComplex</tt>, but you may provide many kernels at once
2003  (using an iterator pair specifying the kernel sequence).
2004  This has the advantage that the forward transform of the input array needs
2005  to be executed only once.
2006  </DL>
2007 
2008  The output arrays must have the same shape as the input arrays. In the "Many" variants of the
2009  convolution functions, the kernels must all have the same shape.
2010 
2011  The origin of the kernel is always assumed to be in the center of the kernel array (precisely,
2012  at the point <tt>floor(kernel.shape() / 2.0)</tt>, except when the half-space format is used, see below).
2013  The function \ref moveDCToUpperLeft() will be called internally to align the kernel with the transformed
2014  input as appropriate.
2015 
2016  If a real input is combined with a real kernel, the kernel is automatically assumed to be defined
2017  in the spatial domain. If a real input is combined with a complex kernel, the kernel is assumed
2018  to be defined in the Fourier domain in half-space format. If the input array is complex, a flag
2019  <tt>fourierDomainKernel</tt> determines where the kernel is defined.
2020 
2021  When the kernel is defined in the spatial domain, the convolution functions will automatically pad
2022  (enlarge) the input array by at least the kernel radius in each direction. The newly added space is
2023  filled according to reflective boundary conditions in order to minimize border artifacts during
2024  convolution. It is thus ensured that convolution in the Fourier domain yields the same results as
2025  convolution in the spatial domain (e.g. when \ref separableConvolveMultiArray() is called with the
2026  same kernel). A little further padding may be added to make sure that the padded array shape
2027  uses integers which have only small prime factors, because FFTW is then able to use the fastest
2028  possible algorithms. Any padding is automatically removed from the result arrays before the function
2029  returns.
2030 
2031  When the kernel is defined in the frequency domain, it must be complex-valued, and its shape determines
2032  the shape of the Fourier representation (i.e. the input is padded according to the shape of the kernel).
2033  If we are going to perform a complex-valued convolution, the kernel must be defined for the entire
2034  frequency domain, and its shape directly determines the size of the FFT.
2035 
2036  In contrast, a frequency domain kernel for a real-valued convolution must have symmetry properties
2037  that allow to drop half of the kernel coefficients, as in the
2038  <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C transform</a>.
2039  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
2040  logical shape of the frequency representation (and thus the size of the padded input). The origin
2041  of the kernel must be at the point
2042  <tt>(0, floor(fourier_shape[0] / 2.0), ..., floor(fourier_shape[N-1] / 2.0))</tt>
2043  (i.e. as in a regular kernel except for the first dimension).
2044 
2045  The <tt>Real</tt> type in the declarations can be <tt>double</tt>, <tt>float</tt>, and
2046  <tt>long double</tt>. Your program must always link against <tt>libfftw3</tt>. If you use
2047  <tt>float</tt> or <tt>long double</tt> arrays, you must <i>additionally</i> link against
2048  <tt>libfftw3f</tt> and <tt>libfftw3l</tt> respectively.
2049 
2050  The Fourier transform functions internally create <a href="http://www.fftw.org/doc/Using-Plans.html">FFTW plans</a>
2051  which control the algorithm details. The plans are creates with the flag <tt>FFTW_ESTIMATE</tt>, i.e.
2052  optimal settings are guessed or read from saved "wisdom" files. If you need more control over planning,
2053  you can use the class \ref FFTWConvolvePlan.
2054 
2055  See also \ref applyFourierFilter() for corresponding functionality on the basis of the
2056  old image iterator interface.
2057 
2058  <b> Declarations:</b>
2059 
2060  Real-valued convolution with kernel in the spatial domain:
2061  \code
2062  namespace vigra {
2063  template <unsigned int N, class Real, class C1, class C2, class C3>
2064  void
2065  convolveFFT(MultiArrayView<N, Real, C1> in,
2066  MultiArrayView<N, Real, C2> kernel,
2067  MultiArrayView<N, Real, C3> out);
2068  }
2069  \endcode
2070 
2071  Real-valued convolution with kernel in the Fourier domain (half-space format):
2072  \code
2073  namespace vigra {
2074  template <unsigned int N, class Real, class C1, class C2, class C3>
2075  void
2076  convolveFFT(MultiArrayView<N, Real, C1> in,
2077  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2078  MultiArrayView<N, Real, C3> out);
2079  }
2080  \endcode
2081 
2082  Series of real-valued convolutions with kernels in the spatial or Fourier domain
2083  (the kernel and out sequences must have the same length):
2084  \code
2085  namespace vigra {
2086  template <unsigned int N, class Real, class C1,
2087  class KernelIterator, class OutIterator>
2088  void
2089  convolveFFTMany(MultiArrayView<N, Real, C1> in,
2090  KernelIterator kernels, KernelIterator kernelsEnd,
2091  OutIterator outs);
2092  }
2093  \endcode
2094 
2095  Complex-valued convolution (parameter <tt>fourierDomainKernel</tt> determines if
2096  the kernel is defined in the spatial or Fourier domain):
2097  \code
2098  namespace vigra {
2099  template <unsigned int N, class Real, class C1, class C2, class C3>
2100  void
2101  convolveFFTComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in,
2102  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2103  MultiArrayView<N, FFTWComplex<Real>, C3> out,
2104  bool fourierDomainKernel);
2105  }
2106  \endcode
2107 
2108  Series of complex-valued convolutions (parameter <tt>fourierDomainKernel</tt>
2109  determines if the kernels are defined in the spatial or Fourier domain,
2110  the kernel and out sequences must have the same length):
2111  \code
2112  namespace vigra {
2113  template <unsigned int N, class Real, class C1,
2114  class KernelIterator, class OutIterator>
2115  void
2116  convolveFFTComplexMany(MultiArrayView<N, FFTWComplex<Real>, C1> in,
2117  KernelIterator kernels, KernelIterator kernelsEnd,
2118  OutIterator outs,
2119  bool fourierDomainKernel);
2120  }
2121  \endcode
2122 
2123  <b> Usage:</b>
2124 
2125  <b>\#include</b> <vigra/multi_fft.hxx><br>
2126  Namespace: vigra
2127 
2128  \code
2129  // convolve real array with a Gaussian (sigma=1) defined in the spatial domain
2130  // (implicitly uses padding by at least 4 pixels)
2131  MultiArray<2, double> src(Shape2(w, h)), dest(Shape2(w,h));
2132 
2133  MultiArray<2, double> spatial_kernel(Shape2(9, 9));
2134  Gaussian<double> gauss(1.0);
2135 
2136  for(int y=0; y<9; ++y)
2137  for(int x=0; x<9; ++x)
2138  spatial_kernel(x, y) = gauss(x-4.0)*gauss(y-4.0);
2139 
2140  convolveFFT(src, spatial_kernel, dest);
2141 
2142  // convolve real array with a Gaussian (sigma=1) defined in the Fourier domain
2143  // (uses no padding, because the kernel size corresponds to the input size)
2144  MultiArray<2, FFTWComplex<double> > fourier_kernel(fftwCorrespondingShapeR2C(src.shape()));
2145  int y0 = h / 2;
2146 
2147  for(int y=0; y<fourier_kernel.shape(1); ++y)
2148  for(int x=0; x<fourier_kernel.shape(0); ++x)
2149  fourier_kernel(x, y) = exp(-0.5*sq(x / double(w))) * exp(-0.5*sq((y-y0)/double(h)));
2150 
2151  convolveFFT(src, fourier_kernel, dest);
2152  \endcode
2153 */
2154 doxygen_overloaded_function(template <...> void convolveFFT)
2155 
2156 template <unsigned int N, class Real, class C1, class C2, class C3>
2157 void
2158 convolveFFT(MultiArrayView<N, Real, C1> in,
2159  MultiArrayView<N, Real, C2> kernel,
2160  MultiArrayView<N, Real, C3> out)
2161 {
2162  FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out);
2163 }
2164 
2165 template <unsigned int N, class Real, class C1, class C2, class C3>
2166 void
2167 convolveFFT(MultiArrayView<N, Real, C1> in,
2168  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2169  MultiArrayView<N, Real, C3> out)
2170 {
2171  FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out);
2172 }
2173 
2174 /** \brief Convolve a complex-valued array by means of the Fourier transform.
2175 
2176  See \ref convolveFFT() for details.
2177 */
2178 doxygen_overloaded_function(template <...> void convolveFFTComplex)
2179 
2180 template <unsigned int N, class Real, class C1, class C2, class C3>
2181 void
2182 convolveFFTComplex(MultiArrayView<N, FFTWComplex<Real>, C1> in,
2183  MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2184  MultiArrayView<N, FFTWComplex<Real>, C3> out,
2185  bool fourierDomainKernel)
2186 {
2187  FFTWConvolvePlan<N, Real>(in, kernel, out, fourierDomainKernel).execute(in, kernel, out);
2188 }
2189 
2190 /** \brief Convolve a real-valued array with a sequence of kernels by means of the Fourier transform.
2191 
2192  See \ref convolveFFT() for details.
2193 */
2194 doxygen_overloaded_function(template <...> void convolveFFTMany)
2195 
2196 template <unsigned int N, class Real, class C1,
2197  class KernelIterator, class OutIterator>
2198 void
2199 convolveFFTMany(MultiArrayView<N, Real, C1> in,
2200  KernelIterator kernels, KernelIterator kernelsEnd,
2201  OutIterator outs)
2202 {
2203  FFTWConvolvePlan<N, Real> plan;
2204  plan.initMany(in, kernels, kernelsEnd, outs);
2205  plan.executeMany(in, kernels, kernelsEnd, outs);
2206 }
2207 
2208 /** \brief Convolve a complex-valued array with a sequence of kernels by means of the Fourier transform.
2209 
2210  See \ref convolveFFT() for details.
2211 */
2212 doxygen_overloaded_function(template <...> void convolveFFTComplexMany)
2213 
2214 template <unsigned int N, class Real, class C1,
2215  class KernelIterator, class OutIterator>
2216 void
2217 convolveFFTComplexMany(MultiArrayView<N, FFTWComplex<Real>, C1> in,
2218  KernelIterator kernels, KernelIterator kernelsEnd,
2219  OutIterator outs,
2220  bool fourierDomainKernel)
2221 {
2222  FFTWConvolvePlan<N, Real> plan;
2223  plan.initMany(in, kernels, kernelsEnd, outs, fourierDomainKernel);
2224  plan.executeMany(in, kernels, kernelsEnd, outs);
2225 }
2226 
2227 //@}
2228 
2229 } // namespace vigra
2230 
2231 #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 (Wed Sep 26 2012)