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

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

html generated using doxygen and Python
vigra 1.11.1