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

multi_fft.hxx VIGRA

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

© 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.10.0 (Fri Feb 21 2014)