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

© 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.11.1 (Fri May 19 2017)