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

multi_fft.hxx VIGRA

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

© 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.0