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

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

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

html generated using doxygen and Python
vigra 1.9.0 (Tue Sep 24 2013)