36 #ifndef VIGRA_MULTI_FFT_HXX
37 #define VIGRA_MULTI_FFT_HXX
40 #include "multi_array.hxx"
41 #include "navigator.hxx"
42 #include "copyimage.hxx"
58 template <
unsigned int N,
class T,
class C>
59 void moveDCToCenterImpl(MultiArrayView<N, T, C> a,
unsigned int startDimension)
61 typedef typename MultiArrayView<N, T, C>::traverser Traverser;
62 typedef MultiArrayNavigator<Traverser, N> Navigator;
63 typedef typename Navigator::iterator Iterator;
65 for(
unsigned int d = startDimension; d < N; ++d)
67 Navigator nav(a.traverser_begin(), a.shape(), d);
69 for( ; nav.hasMore(); nav++ )
71 Iterator i = nav.begin();
72 int s = nav.end() - i;
77 for(
int k=0; k<s2; ++k)
79 std::swap(i[k], i[k+s2]);
85 for(
int k=0; k<s2; ++k)
96 template <
unsigned int N,
class T,
class C>
97 void moveDCToUpperLeftImpl(MultiArrayView<N, T, C> a,
unsigned int startDimension)
99 typedef typename MultiArrayView<N, T, C>::traverser Traverser;
100 typedef MultiArrayNavigator<Traverser, N> Navigator;
101 typedef typename Navigator::iterator Iterator;
103 for(
unsigned int d = startDimension; d < N; ++d)
105 Navigator nav(a.traverser_begin(), a.shape(), d);
107 for( ; nav.hasMore(); nav++ )
109 Iterator i = nav.begin();
110 int s = nav.end() - i;
115 for(
int k=0; k<s2; ++k)
117 std::swap(i[k], i[k+s2]);
123 for(
int k=s2; k>0; --k)
142 template <
unsigned int N,
class T,
class C>
145 detail::moveDCToCenterImpl(a, 0);
148 template <
unsigned int N,
class T,
class C>
151 detail::moveDCToUpperLeftImpl(a, 0);
154 template <
unsigned int N,
class T,
class C>
155 inline void moveDCToHalfspaceCenter(MultiArrayView<N, T, C> a)
157 detail::moveDCToCenterImpl(a, 1);
160 template <
unsigned int N,
class T,
class C>
161 inline void moveDCToHalfspaceUpperLeft(MultiArrayView<N, T, C> a)
163 detail::moveDCToUpperLeftImpl(a, 1);
170 fftwPlanCreate(
unsigned int N,
int* shape,
171 FFTWComplex<double> * in,
int* instrides,
int instep,
172 FFTWComplex<double> * out,
int* outstrides,
int outstep,
173 int sign,
unsigned int planner_flags)
175 return fftw_plan_many_dft(N, shape, 1,
176 (fftw_complex *)in, instrides, instep, 0,
177 (fftw_complex *)out, outstrides, outstep, 0,
178 sign, planner_flags);
182 fftwPlanCreate(
unsigned int N,
int* shape,
183 double * in,
int* instrides,
int instep,
184 FFTWComplex<double> * out,
int* outstrides,
int outstep,
185 int ,
unsigned int planner_flags)
187 return fftw_plan_many_dft_r2c(N, shape, 1,
188 in, instrides, instep, 0,
189 (fftw_complex *)out, outstrides, outstep, 0,
194 fftwPlanCreate(
unsigned int N,
int* shape,
195 FFTWComplex<double> * in,
int* instrides,
int instep,
196 double * out,
int* outstrides,
int outstep,
197 int ,
unsigned int planner_flags)
199 return fftw_plan_many_dft_c2r(N, shape, 1,
200 (fftw_complex *)in, instrides, instep, 0,
201 out, outstrides, outstep, 0,
206 fftwPlanCreate(
unsigned int N,
int* shape,
207 FFTWComplex<float> * in,
int* instrides,
int instep,
208 FFTWComplex<float> * out,
int* outstrides,
int outstep,
209 int sign,
unsigned int planner_flags)
211 return fftwf_plan_many_dft(N, shape, 1,
212 (fftwf_complex *)in, instrides, instep, 0,
213 (fftwf_complex *)out, outstrides, outstep, 0,
214 sign, planner_flags);
218 fftwPlanCreate(
unsigned int N,
int* shape,
219 float * in,
int* instrides,
int instep,
220 FFTWComplex<float> * out,
int* outstrides,
int outstep,
221 int ,
unsigned int planner_flags)
223 return fftwf_plan_many_dft_r2c(N, shape, 1,
224 in, instrides, instep, 0,
225 (fftwf_complex *)out, outstrides, outstep, 0,
230 fftwPlanCreate(
unsigned int N,
int* shape,
231 FFTWComplex<float> * in,
int* instrides,
int instep,
232 float * out,
int* outstrides,
int outstep,
233 int ,
unsigned int planner_flags)
235 return fftwf_plan_many_dft_c2r(N, shape, 1,
236 (fftwf_complex *)in, instrides, instep, 0,
237 out, outstrides, outstep, 0,
242 fftwPlanCreate(
unsigned int N,
int* shape,
243 FFTWComplex<long double> * in,
int* instrides,
int instep,
244 FFTWComplex<long double> * out,
int* outstrides,
int outstep,
245 int sign,
unsigned int planner_flags)
247 return fftwl_plan_many_dft(N, shape, 1,
248 (fftwl_complex *)in, instrides, instep, 0,
249 (fftwl_complex *)out, outstrides, outstep, 0,
250 sign, planner_flags);
254 fftwPlanCreate(
unsigned int N,
int* shape,
255 long double * in,
int* instrides,
int instep,
256 FFTWComplex<long double> * out,
int* outstrides,
int outstep,
257 int ,
unsigned int planner_flags)
259 return fftwl_plan_many_dft_r2c(N, shape, 1,
260 in, instrides, instep, 0,
261 (fftwl_complex *)out, outstrides, outstep, 0,
266 fftwPlanCreate(
unsigned int N,
int* shape,
267 FFTWComplex<long double> * in,
int* instrides,
int instep,
268 long double * out,
int* outstrides,
int outstep,
269 int ,
unsigned int planner_flags)
271 return fftwl_plan_many_dft_c2r(N, shape, 1,
272 (fftwl_complex *)in, instrides, instep, 0,
273 out, outstrides, outstep, 0,
277 inline void fftwPlanDestroy(fftw_plan plan)
280 fftw_destroy_plan(plan);
283 inline void fftwPlanDestroy(fftwf_plan plan)
286 fftwf_destroy_plan(plan);
289 inline void fftwPlanDestroy(fftwl_plan plan)
292 fftwl_destroy_plan(plan);
296 fftwPlanExecute(fftw_plan plan)
302 fftwPlanExecute(fftwf_plan plan)
308 fftwPlanExecute(fftwl_plan plan)
314 fftwPlanExecute(fftw_plan plan, FFTWComplex<double> * in, FFTWComplex<double> * out)
316 fftw_execute_dft(plan, (fftw_complex *)in, (fftw_complex *)out);
320 fftwPlanExecute(fftw_plan plan,
double * in, FFTWComplex<double> * out)
322 fftw_execute_dft_r2c(plan, in, (fftw_complex *)out);
326 fftwPlanExecute(fftw_plan plan, FFTWComplex<double> * in,
double * out)
328 fftw_execute_dft_c2r(plan, (fftw_complex *)in, out);
332 fftwPlanExecute(fftwf_plan plan, FFTWComplex<float> * in, FFTWComplex<float> * out)
334 fftwf_execute_dft(plan, (fftwf_complex *)in, (fftwf_complex *)out);
338 fftwPlanExecute(fftwf_plan plan,
float * in, FFTWComplex<float> * out)
340 fftwf_execute_dft_r2c(plan, in, (fftwf_complex *)out);
344 fftwPlanExecute(fftwf_plan plan, FFTWComplex<float> * in,
float * out)
346 fftwf_execute_dft_c2r(plan, (fftwf_complex *)in, out);
350 fftwPlanExecute(fftwl_plan plan, FFTWComplex<long double> * in, FFTWComplex<long double> * out)
352 fftwl_execute_dft(plan, (fftwl_complex *)in, (fftwl_complex *)out);
356 fftwPlanExecute(fftwl_plan plan,
long double * in, FFTWComplex<long double> * out)
358 fftwl_execute_dft_r2c(plan, in, (fftwl_complex *)out);
362 fftwPlanExecute(fftwl_plan plan, FFTWComplex<long double> * in,
long double * out)
364 fftwl_execute_dft_c2r(plan, (fftwl_complex *)in, out);
368 int fftwPaddingSize(
int s)
374 static const int size = 1330;
375 static int goodSizes[size] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
376 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 35, 36, 39, 40, 42, 44, 45, 48,
377 49, 50, 52, 54, 55, 56, 60, 63, 64, 65, 66, 70, 72, 75, 77, 78, 80, 81,
378 84, 88, 90, 91, 96, 98, 99, 100, 104, 105, 108, 110, 112, 117, 120, 125,
379 126, 128, 130, 132, 135, 140, 144, 147, 150, 154, 156, 160, 162, 165,
380 168, 175, 176, 180, 182, 189, 192, 195, 196, 198, 200, 208, 210, 216,
381 220, 224, 225, 231, 234, 240, 243, 245, 250, 252, 256, 260, 264, 270,
382 273, 275, 280, 288, 294, 297, 300, 308, 312, 315, 320, 324, 325, 330,
383 336, 343, 350, 351, 352, 360, 364, 375, 378, 384, 385, 390, 392, 396,
384 400, 405, 416, 420, 432, 440, 441, 448, 450, 455, 462, 468, 480, 486,
385 490, 495, 500, 504, 512, 520, 525, 528, 539, 540, 546, 550, 560, 567,
386 576, 585, 588, 594, 600, 616, 624, 625, 630, 637, 640, 648, 650, 660,
387 672, 675, 686, 693, 700, 702, 704, 720, 728, 729, 735, 750, 756, 768,
388 770, 780, 784, 792, 800, 810, 819, 825, 832, 840, 864, 875, 880, 882,
389 891, 896, 900, 910, 924, 936, 945, 960, 972, 975, 980, 990, 1000, 1008,
390 1024, 1029, 1040, 1050, 1053, 1056, 1078, 1080, 1092, 1100, 1120, 1125,
391 1134, 1152, 1155, 1170, 1176, 1188, 1200, 1215, 1225, 1232, 1248, 1250,
392 1260, 1274, 1280, 1296, 1300, 1320, 1323, 1344, 1350, 1365, 1372, 1375,
393 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470, 1485, 1500, 1512, 1536,
394 1540, 1560, 1568, 1575, 1584, 1600, 1617, 1620, 1625, 1638, 1650, 1664,
395 1680, 1701, 1715, 1728, 1750, 1755, 1760, 1764, 1782, 1792, 1800, 1820,
396 1848, 1872, 1875, 1890, 1911, 1920, 1925, 1944, 1950, 1960, 1980, 2000,
397 2016, 2025, 2048, 2058, 2079, 2080, 2100, 2106, 2112, 2156, 2160, 2184,
398 2187, 2200, 2205, 2240, 2250, 2268, 2275, 2304, 2310, 2340, 2352, 2376,
399 2400, 2401, 2430, 2450, 2457, 2464, 2475, 2496, 2500, 2520, 2548, 2560,
400 2592, 2600, 2625, 2640, 2646, 2673, 2688, 2695, 2700, 2730, 2744, 2750,
401 2772, 2800, 2808, 2816, 2835, 2880, 2912, 2916, 2925, 2940, 2970, 3000,
402 3024, 3072, 3080, 3087, 3120, 3125, 3136, 3150, 3159, 3168, 3185, 3200,
403 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3375, 3402, 3430, 3456, 3465,
404 3500, 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3645, 3675, 3696, 3744,
405 3750, 3773, 3780, 3822, 3840, 3850, 3888, 3900, 3920, 3960, 3969, 4000,
406 4032, 4050, 4095, 4096, 4116, 4125, 4158, 4160, 4200, 4212, 4224, 4312,
407 4320, 4368, 4374, 4375, 4400, 4410, 4455, 4459, 4480, 4500, 4536, 4550,
408 4608, 4620, 4680, 4704, 4725, 4752, 4800, 4802, 4851, 4860, 4875, 4900,
409 4914, 4928, 4950, 4992, 5000, 5040, 5096, 5103, 5120, 5145, 5184, 5200,
410 5250, 5265, 5280, 5292, 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544,
411 5600, 5616, 5625, 5632, 5670, 5733, 5760, 5775, 5824, 5832, 5850, 5880,
412 5940, 6000, 6048, 6075, 6125, 6144, 6160, 6174, 6237, 6240, 6250, 6272,
413 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6561, 6600, 6615,
414 6656, 6720, 6750, 6804, 6825, 6860, 6875, 6912, 6930, 7000, 7020, 7040,
415 7056, 7128, 7168, 7200, 7203, 7280, 7290, 7350, 7371, 7392, 7425, 7488,
416 7500, 7546, 7560, 7644, 7680, 7700, 7776, 7800, 7840, 7875, 7920, 7938,
417 8000, 8019, 8064, 8085, 8100, 8125, 8190, 8192, 8232, 8250, 8316, 8320,
418 8400, 8424, 8448, 8505, 8575, 8624, 8640, 8736, 8748, 8750, 8775, 8800,
419 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9261, 9360, 9375,
420 9408, 9450, 9477, 9504, 9555, 9600, 9604, 9625, 9702, 9720, 9750, 9800,
421 9828, 9856, 9900, 9984, 10000, 10080, 10125, 10192, 10206, 10240, 10290,
422 10368, 10395, 10400, 10500, 10530, 10560, 10584, 10692, 10752, 10780,
423 10800, 10920, 10935, 10976, 11000, 11025, 11088, 11200, 11232, 11250,
424 11264, 11319, 11340, 11375, 11466, 11520, 11550, 11648, 11664, 11700,
425 11760, 11880, 11907, 12000, 12005, 12096, 12150, 12250, 12285, 12288,
426 12320, 12348, 12375, 12474, 12480, 12500, 12544, 12600, 12636, 12672,
427 12740, 12800, 12936, 12960, 13000, 13104, 13122, 13125, 13200, 13230,
428 13312, 13365, 13377, 13440, 13475, 13500, 13608, 13650, 13720, 13750,
429 13824, 13860, 14000, 14040, 14080, 14112, 14175, 14256, 14336, 14400,
430 14406, 14553, 14560, 14580, 14625, 14700, 14742, 14784, 14850, 14976,
431 15000, 15092, 15120, 15288, 15309, 15360, 15400, 15435, 15552, 15600,
432 15625, 15680, 15750, 15795, 15840, 15876, 15925, 16000, 16038, 16128,
433 16170, 16200, 16250, 16380, 16384, 16464, 16500, 16632, 16640, 16800,
434 16807, 16848, 16875, 16896, 17010, 17150, 17199, 17248, 17280, 17325,
435 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000,
436 18144, 18200, 18225, 18375, 18432, 18480, 18522, 18711, 18720, 18750,
437 18816, 18865, 18900, 18954, 19008, 19110, 19200, 19208, 19250, 19404,
438 19440, 19500, 19600, 19656, 19683, 19712, 19800, 19845, 19968, 20000,
439 20160, 20250, 20384, 20412, 20475, 20480, 20580, 20625, 20736, 20790,
440 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560, 21600, 21609,
441 21840, 21870, 21875, 21952, 22000, 22050, 22113, 22176, 22275, 22295,
442 22400, 22464, 22500, 22528, 22638, 22680, 22750, 22932, 23040, 23100,
443 23296, 23328, 23400, 23520, 23625, 23760, 23814, 24000, 24010, 24057,
444 24192, 24255, 24300, 24375, 24500, 24570, 24576, 24640, 24696, 24750,
445 24948, 24960, 25000, 25088, 25200, 25272, 25344, 25480, 25515, 25600,
446 25725, 25872, 25920, 26000, 26208, 26244, 26250, 26325, 26400, 26411,
447 26460, 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440,
448 27500, 27648, 27720, 27783, 28000, 28080, 28125, 28160, 28224, 28350,
449 28431, 28512, 28665, 28672, 28800, 28812, 28875, 29106, 29120, 29160,
450 29250, 29400, 29484, 29568, 29700, 29952, 30000, 30184, 30240, 30375,
451 30576, 30618, 30625, 30720, 30800, 30870, 31104, 31185, 31200, 31213,
452 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256,
453 32340, 32400, 32500, 32760, 32768, 32805, 32928, 33000, 33075, 33264,
454 33280, 33600, 33614, 33696, 33750, 33792, 33957, 34020, 34125, 34300,
455 34375, 34398, 34496, 34560, 34650, 34944, 34992, 35000, 35100, 35200,
456 35280, 35640, 35672, 35721, 35840, 36000, 36015, 36288, 36400, 36450,
457 36750, 36855, 36864, 36960, 37044, 37125, 37422, 37440, 37500, 37632,
458 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808, 38880,
459 39000, 39200, 39312, 39366, 39375, 39424, 39600, 39690, 39936, 40000,
460 40095, 40131, 40320, 40425, 40500, 40625, 40768, 40824, 40950, 40960,
461 41160, 41250, 41472, 41580, 41600, 42000, 42120, 42240, 42336, 42525,
462 42768, 42875, 43008, 43120, 43200, 43218, 43659, 43680, 43740, 43750,
463 43875, 43904, 44000, 44100, 44226, 44352, 44550, 44590, 44800, 44928,
464 45000, 45056, 45276, 45360, 45500, 45864, 45927, 46080, 46200, 46305,
465 46592, 46656, 46800, 46875, 47040, 47250, 47385, 47520, 47628, 47775,
466 48000, 48020, 48114, 48125, 48384, 48510, 48600, 48750, 49000, 49140,
467 49152, 49280, 49392, 49500, 49896, 49920, 50000, 50176, 50400, 50421,
468 50544, 50625, 50688, 50960, 51030, 51200, 51450, 51597, 51744, 51840,
469 51975, 52000, 52416, 52488, 52500, 52650, 52800, 52822, 52920, 53248,
470 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54675, 54880, 55000,
471 55125, 55296, 55440, 55566, 56000, 56133, 56160, 56250, 56320, 56448,
472 56595, 56700, 56862, 56875, 57024, 57330, 57344, 57600, 57624, 57750,
473 58212, 58240, 58320, 58500, 58800, 58968, 59049, 59136, 59400, 59535,
474 59904, 60000, 60025, 60368, 60480, 60750, 61152, 61236, 61250, 61425,
475 61440, 61600, 61740, 61875, 62208, 62370, 62400, 62426, 62500, 62720,
476 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512, 64680, 64800,
477 64827, 65000, 65520, 65536, 65610, 65625, 65856, 66000, 66150, 66339,
478 66528, 66560, 66825, 66885, 67200, 67228, 67375, 67392, 67500, 67584,
479 67914, 68040, 68250, 68600, 68750, 68796, 68992, 69120, 69300, 69888,
480 69984, 70000, 70200, 70400, 70560, 70875, 71280, 71344, 71442, 71680,
481 72000, 72030, 72171, 72576, 72765, 72800, 72900, 73125, 73500, 73710,
482 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264, 75460, 75600,
483 75816, 76032, 76440, 76545, 76800, 76832, 77000, 77175, 77616, 77760,
484 78000, 78125, 78400, 78624, 78732, 78750, 78848, 78975, 79200, 79233,
485 79380, 79625, 79872, 80000, 80190, 80262, 80640, 80850, 81000, 81250,
486 81536, 81648, 81900, 81920, 82320, 82500, 82944, 83160, 83200, 83349,
487 84000, 84035, 84240, 84375, 84480, 84672, 85050, 85293, 85536, 85750,
488 85995, 86016, 86240, 86400, 86436, 86625, 87318, 87360, 87480, 87500,
489 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856,
490 90000, 90112, 90552, 90720, 91000, 91125, 91728, 91854, 91875, 92160,
491 92400, 92610, 93184, 93312, 93555, 93600, 93639, 93750, 94080, 94325,
492 94500, 94770, 95040, 95256, 95550, 96000, 96040, 96228, 96250, 96768,
493 97020, 97200, 97500, 98000, 98280, 98304, 98415, 98560, 98784, 99000,
494 99225, 99792, 99840 };
496 if(s <= 0 || s > goodSizes[size-1])
498 return *std::upper_bound(goodSizes, goodSizes+size, s, std::less_equal<int>());
502 int fftwEvenPaddingSize(
int s)
508 static const int size = 1063;
509 static int goodSizes[size] = { 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22,
510 24, 26, 28, 30, 32, 36, 40, 42, 44, 48, 50, 52, 54, 56, 60, 64, 66, 70,
511 72, 78, 80, 84, 88, 90, 96, 98, 100, 104, 108, 110, 112, 120, 126, 128,
512 130, 132, 140, 144, 150, 154, 156, 160, 162, 168, 176, 180, 182, 192,
513 196, 198, 200, 208, 210, 216, 220, 224, 234, 240, 250, 252, 256, 260,
514 264, 270, 280, 288, 294, 300, 308, 312, 320, 324, 330, 336, 350, 352,
515 360, 364, 378, 384, 390, 392, 396, 400, 416, 420, 432, 440, 448, 450,
516 462, 468, 480, 486, 490, 500, 504, 512, 520, 528, 540, 546, 550, 560,
517 576, 588, 594, 600, 616, 624, 630, 640, 648, 650, 660, 672, 686, 700,
518 702, 704, 720, 728, 750, 756, 768, 770, 780, 784, 792, 800, 810, 832,
519 840, 864, 880, 882, 896, 900, 910, 924, 936, 960, 972, 980, 990, 1000,
520 1008, 1024, 1040, 1050, 1056, 1078, 1080, 1092, 1100, 1120, 1134, 1152,
521 1170, 1176, 1188, 1200, 1232, 1248, 1250, 1260, 1274, 1280, 1296, 1300,
522 1320, 1344, 1350, 1372, 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470,
523 1500, 1512, 1536, 1540, 1560, 1568, 1584, 1600, 1620, 1638, 1650, 1664,
524 1680, 1728, 1750, 1760, 1764, 1782, 1792, 1800, 1820, 1848, 1872, 1890,
525 1920, 1944, 1950, 1960, 1980, 2000, 2016, 2048, 2058, 2080, 2100, 2106,
526 2112, 2156, 2160, 2184, 2200, 2240, 2250, 2268, 2304, 2310, 2340, 2352,
527 2376, 2400, 2430, 2450, 2464, 2496, 2500, 2520, 2548, 2560, 2592, 2600,
528 2640, 2646, 2688, 2700, 2730, 2744, 2750, 2772, 2800, 2808, 2816, 2880,
529 2912, 2916, 2940, 2970, 3000, 3024, 3072, 3080, 3120, 3136, 3150, 3168,
530 3200, 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3402, 3430, 3456, 3500,
531 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3696, 3744, 3750, 3780, 3822,
532 3840, 3850, 3888, 3900, 3920, 3960, 4000, 4032, 4050, 4096, 4116, 4158,
533 4160, 4200, 4212, 4224, 4312, 4320, 4368, 4374, 4400, 4410, 4480, 4500,
534 4536, 4550, 4608, 4620, 4680, 4704, 4752, 4800, 4802, 4860, 4900, 4914,
535 4928, 4950, 4992, 5000, 5040, 5096, 5120, 5184, 5200, 5250, 5280, 5292,
536 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544, 5600, 5616, 5632, 5670,
537 5760, 5824, 5832, 5850, 5880, 5940, 6000, 6048, 6144, 6160, 6174, 6240,
538 6250, 6272, 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6600,
539 6656, 6720, 6750, 6804, 6860, 6912, 6930, 7000, 7020, 7040, 7056, 7128,
540 7168, 7200, 7280, 7290, 7350, 7392, 7488, 7500, 7546, 7560, 7644, 7680,
541 7700, 7776, 7800, 7840, 7920, 7938, 8000, 8064, 8100, 8190, 8192, 8232,
542 8250, 8316, 8320, 8400, 8424, 8448, 8624, 8640, 8736, 8748, 8750, 8800,
543 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9360, 9408, 9450,
544 9504, 9600, 9604, 9702, 9720, 9750, 9800, 9828, 9856, 9900, 9984, 10000,
545 10080, 10192, 10206, 10240, 10290, 10368, 10400, 10500, 10530, 10560,
546 10584, 10692, 10752, 10780, 10800, 10920, 10976, 11000, 11088, 11200,
547 11232, 11250, 11264, 11340, 11466, 11520, 11550, 11648, 11664, 11700,
548 11760, 11880, 12000, 12096, 12150, 12250, 12288, 12320, 12348, 12474,
549 12480, 12500, 12544, 12600, 12636, 12672, 12740, 12800, 12936, 12960,
550 13000, 13104, 13122, 13200, 13230, 13312, 13440, 13500, 13608, 13650,
551 13720, 13750, 13824, 13860, 14000, 14040, 14080, 14112, 14256, 14336,
552 14400, 14406, 14560, 14580, 14700, 14742, 14784, 14850, 14976, 15000,
553 15092, 15120, 15288, 15360, 15400, 15552, 15600, 15680, 15750, 15840,
554 15876, 16000, 16038, 16128, 16170, 16200, 16250, 16380, 16384, 16464,
555 16500, 16632, 16640, 16800, 16848, 16896, 17010, 17150, 17248, 17280,
556 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000,
557 18144, 18200, 18432, 18480, 18522, 18720, 18750, 18816, 18900, 18954,
558 19008, 19110, 19200, 19208, 19250, 19404, 19440, 19500, 19600, 19656,
559 19712, 19800, 19968, 20000, 20160, 20250, 20384, 20412, 20480, 20580,
560 20736, 20790, 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560,
561 21600, 21840, 21870, 21952, 22000, 22050, 22176, 22400, 22464, 22500,
562 22528, 22638, 22680, 22750, 22932, 23040, 23100, 23296, 23328, 23400,
563 23520, 23760, 23814, 24000, 24010, 24192, 24300, 24500, 24570, 24576,
564 24640, 24696, 24750, 24948, 24960, 25000, 25088, 25200, 25272, 25344,
565 25480, 25600, 25872, 25920, 26000, 26208, 26244, 26250, 26400, 26460,
566 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440, 27500,
567 27648, 27720, 28000, 28080, 28160, 28224, 28350, 28512, 28672, 28800,
568 28812, 29106, 29120, 29160, 29250, 29400, 29484, 29568, 29700, 29952,
569 30000, 30184, 30240, 30576, 30618, 30720, 30800, 30870, 31104, 31200,
570 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256,
571 32340, 32400, 32500, 32760, 32768, 32928, 33000, 33264, 33280, 33600,
572 33614, 33696, 33750, 33792, 34020, 34300, 34398, 34496, 34560, 34650,
573 34944, 34992, 35000, 35100, 35200, 35280, 35640, 35672, 35840, 36000,
574 36288, 36400, 36450, 36750, 36864, 36960, 37044, 37422, 37440, 37500,
575 37632, 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808,
576 38880, 39000, 39200, 39312, 39366, 39424, 39600, 39690, 39936, 40000,
577 40320, 40500, 40768, 40824, 40950, 40960, 41160, 41250, 41472, 41580,
578 41600, 42000, 42120, 42240, 42336, 42768, 43008, 43120, 43200, 43218,
579 43680, 43740, 43750, 43904, 44000, 44100, 44226, 44352, 44550, 44590,
580 44800, 44928, 45000, 45056, 45276, 45360, 45500, 45864, 46080, 46200,
581 46592, 46656, 46800, 47040, 47250, 47520, 47628, 48000, 48020, 48114,
582 48384, 48510, 48600, 48750, 49000, 49140, 49152, 49280, 49392, 49500,
583 49896, 49920, 50000, 50176, 50400, 50544, 50688, 50960, 51030, 51200,
584 51450, 51744, 51840, 52000, 52416, 52488, 52500, 52650, 52800, 52822,
585 52920, 53248, 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54880,
586 55000, 55296, 55440, 55566, 56000, 56160, 56250, 56320, 56448, 56700,
587 56862, 57024, 57330, 57344, 57600, 57624, 57750, 58212, 58240, 58320,
588 58500, 58800, 58968, 59136, 59400, 59904, 60000, 60368, 60480, 60750,
589 61152, 61236, 61250, 61440, 61600, 61740, 62208, 62370, 62400, 62426,
590 62500, 62720, 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512,
591 64680, 64800, 65000, 65520, 65536, 65610, 65856, 66000, 66150, 66528,
592 66560, 67200, 67228, 67392, 67500, 67584, 67914, 68040, 68250, 68600,
593 68750, 68796, 68992, 69120, 69300, 69888, 69984, 70000, 70200, 70400,
594 70560, 71280, 71344, 71442, 71680, 72000, 72030, 72576, 72800, 72900,
595 73500, 73710, 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264,
596 75460, 75600, 75816, 76032, 76440, 76800, 76832, 77000, 77616, 77760,
597 78000, 78400, 78624, 78732, 78750, 78848, 79200, 79380, 79872, 80000,
598 80190, 80262, 80640, 80850, 81000, 81250, 81536, 81648, 81900, 81920,
599 82320, 82500, 82944, 83160, 83200, 84000, 84240, 84480, 84672, 85050,
600 85536, 85750, 86016, 86240, 86400, 86436, 87318, 87360, 87480, 87500,
601 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856,
602 90000, 90112, 90552, 90720, 91000, 91728, 91854, 92160, 92400, 92610,
603 93184, 93312, 93600, 93750, 94080, 94500, 94770, 95040, 95256, 95550,
604 96000, 96040, 96228, 96250, 96768, 97020, 97200, 97500, 98000, 98280,
605 98304, 98560, 98784, 99000, 99792, 99840 };
607 if(s <= 0 || s > goodSizes[size-1])
609 return *std::upper_bound(goodSizes, goodSizes+size, s, std::less_equal<int>());
613 struct FFTEmbedKernel
615 template <
unsigned int N,
class Real,
class C,
class Shape>
617 exec(MultiArrayView<N, Real, C> & out, Shape
const & kernelShape,
618 Shape & srcPoint, Shape & destPoint,
bool copyIt)
620 for(srcPoint[M]=0; srcPoint[M]<kernelShape[M]; ++srcPoint[M])
622 if(srcPoint[M] < (kernelShape[M] + 1) / 2)
624 destPoint[M] = srcPoint[M];
628 destPoint[M] = srcPoint[M] + out.shape(M) - kernelShape[M];
631 FFTEmbedKernel<M-1>::exec(out, kernelShape, srcPoint, destPoint, copyIt);
637 struct FFTEmbedKernel<0>
639 template <
unsigned int N,
class Real,
class C,
class Shape>
641 exec(MultiArrayView<N, Real, C> & out, Shape
const & kernelShape,
642 Shape & srcPoint, Shape & destPoint,
bool copyIt)
644 for(srcPoint[0]=0; srcPoint[0]<kernelShape[0]; ++srcPoint[0])
646 if(srcPoint[0] < (kernelShape[0] + 1) / 2)
648 destPoint[0] = srcPoint[0];
652 destPoint[0] = srcPoint[0] + out.shape(0) - kernelShape[0];
657 out[destPoint] = out[srcPoint];
664 template <
unsigned int N,
class Real,
class C1,
class C2>
666 fftEmbedKernel(MultiArrayView<N, Real, C1> kernel,
667 MultiArrayView<N, Real, C2> out,
670 typedef typename MultiArrayShape<N>::type Shape;
672 MultiArrayView<N, Real, C2> kout = out.subarray(Shape(), kernel.shape());
680 Shape srcPoint, destPoint;
681 FFTEmbedKernel<(int)N-1>::exec(out, kernel.shape(), srcPoint, destPoint,
false);
684 template <
unsigned int N,
class Real,
class C1,
class C2>
686 fftEmbedArray(MultiArrayView<N, Real, C1> in,
687 MultiArrayView<N, Real, C2> out)
689 typedef typename MultiArrayShape<N>::type Shape;
691 Shape diff = out.shape() - in.shape(),
693 rightDiff = diff - leftDiff,
694 right = in.shape() + leftDiff;
696 out.subarray(leftDiff, right) = in;
698 typedef typename MultiArrayView<N, Real, C2>::traverser Traverser;
699 typedef MultiArrayNavigator<Traverser, N> Navigator;
700 typedef typename Navigator::iterator Iterator;
702 for(
unsigned int d = 0; d < N; ++d)
704 Navigator nav(out.traverser_begin(), out.shape(), d);
706 for( ; nav.hasMore(); nav++ )
708 Iterator i = nav.begin();
709 for(
int k=1; k<=leftDiff[d]; ++k)
710 i[leftDiff[d] - k] = i[leftDiff[d] + k];
711 for(
int k=0; k<rightDiff[d]; ++k)
712 i[right[d] + k] = i[right[d] - k - 2];
719 template <
class T,
int N>
721 fftwBestPaddedShape(TinyVector<T, N> shape)
723 for(
unsigned int k=0; k<N; ++k)
724 shape[k] = detail::fftwPaddingSize(shape[k]);
728 template <
class T,
int N>
730 fftwBestPaddedShapeR2C(TinyVector<T, N> shape)
732 shape[0] = detail::fftwEvenPaddingSize(shape[0]);
733 for(
unsigned int k=1; k<N; ++k)
734 shape[k] = detail::fftwPaddingSize(shape[k]);
748 template <
class T,
int N>
752 shape[0] = shape[0] / 2 + 1;
756 template <
class T,
int N>
758 fftwCorrespondingShapeC2R(TinyVector<T, N> shape,
bool oddDimension0 =
false)
760 shape[0] = oddDimension0
761 ? (shape[0] - 1) * 2 + 1
762 : (shape[0] - 1) * 2;
799 template <
unsigned int N,
class Real =
double>
803 typedef typename FFTWReal2Complex<Real>::plan_type PlanType;
807 Shape shape, instrides, outstrides;
827 template <
class C1,
class C2>
830 int SIGN,
unsigned int planner_flags = FFTW_ESTIMATE)
833 init(in, out, SIGN, planner_flags);
847 template <
class C1,
class C2>
850 unsigned int planner_flags = FFTW_ESTIMATE)
853 init(in, out, planner_flags);
867 template <
class C1,
class C2>
870 unsigned int planner_flags = FFTW_ESTIMATE)
873 init(in, out, planner_flags);
884 instrides.swap(o.instrides);
885 outstrides.swap(o.outstrides);
898 instrides.swap(o.instrides);
899 outstrides.swap(o.outstrides);
910 detail::fftwPlanDestroy(plan);
917 template <
class C1,
class C2>
920 int SIGN,
unsigned int planner_flags = FFTW_ESTIMATE)
922 vigra_precondition(in.strideOrdering() == out.strideOrdering(),
923 "FFTWPlan.init(): input and output must have the same stride ordering.");
925 initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(),
926 SIGN, planner_flags);
933 template <
class C1,
class C2>
936 unsigned int planner_flags = FFTW_ESTIMATE)
939 "FFTWPlan.init(): input and output must have the same stride ordering.");
942 FFTW_FORWARD, planner_flags);
949 template <
class C1,
class C2>
952 unsigned int planner_flags = FFTW_ESTIMATE)
955 "FFTWPlan.init(): input and output must have the same stride ordering.");
958 FFTW_BACKWARD, planner_flags);
968 template <
class C1,
class C2>
972 executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending());
982 template <
class C1,
class C2>
996 template <
class C1,
class C2>
1005 template <
class MI,
class MO>
1006 void initImpl(MI ins, MO outs,
int SIGN,
unsigned int planner_flags);
1008 template <
class MI,
class MO>
1009 void executeImpl(MI ins, MO outs)
const;
1014 vigra_precondition(in.shape() == out.shape(),
1015 "FFTWPlan.init(): input and output must have the same shape.");
1018 void checkShapes(MultiArrayView<N, Real, StridedArrayTag> ins,
1019 MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> outs)
const
1021 for(
int k=0; k<(int)N-1; ++k)
1022 vigra_precondition(ins.shape(k) == outs.shape(k),
1023 "FFTWPlan.init(): input and output must have matching shapes.");
1024 vigra_precondition(ins.shape(N-1) / 2 + 1 == outs.shape(N-1),
1025 "FFTWPlan.init(): input and output must have matching shapes.");
1028 void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> ins,
1029 MultiArrayView<N, Real, StridedArrayTag> outs)
const
1031 for(
int k=0; k<(int)N-1; ++k)
1032 vigra_precondition(ins.shape(k) == outs.shape(k),
1033 "FFTWPlan.init(): input and output must have matching shapes.");
1034 vigra_precondition(outs.shape(N-1) / 2 + 1 == ins.shape(N-1),
1035 "FFTWPlan.init(): input and output must have matching shapes.");
1039 template <
unsigned int N,
class Real>
1040 template <
class MI,
class MO>
1042 FFTWPlan<N, Real>::initImpl(MI ins, MO outs,
int SIGN,
unsigned int planner_flags)
1044 checkShapes(ins, outs);
1046 typename MultiArrayShape<N>::type logicalShape(SIGN == FFTW_FORWARD
1050 Shape newShape(logicalShape.begin(), logicalShape.end()),
1051 newIStrides(ins.stride().begin(), ins.stride().end()),
1052 newOStrides(outs.stride().begin(), outs.stride().end()),
1053 itotal(ins.shape().begin(), ins.shape().end()),
1054 ototal(outs.shape().begin(), outs.shape().end());
1056 for(
unsigned int j=1; j<N; ++j)
1058 itotal[j] = ins.stride(j-1) / ins.stride(j);
1059 ototal[j] = outs.stride(j-1) / outs.stride(j);
1062 PlanType newPlan = detail::fftwPlanCreate(N, newShape.begin(),
1063 ins.data(), itotal.begin(), ins.stride(N-1),
1064 outs.data(), ototal.begin(), outs.stride(N-1),
1065 SIGN, planner_flags);
1066 detail::fftwPlanDestroy(plan);
1068 shape.swap(newShape);
1069 instrides.swap(newIStrides);
1070 outstrides.swap(newOStrides);
1074 template <
unsigned int N,
class Real>
1075 template <
class MI,
class MO>
1076 void FFTWPlan<N, Real>::executeImpl(MI ins, MO outs)
const
1078 vigra_precondition(plan != 0,
"FFTWPlan::execute(): plan is NULL.");
1080 typename MultiArrayShape<N>::type lshape(sign == FFTW_FORWARD
1084 vigra_precondition((lshape == TinyVectorView<int, N>(shape.data())),
1085 "FFTWPlan::execute(): shape mismatch between plan and data.");
1086 vigra_precondition((ins.stride() == TinyVectorView<int, N>(instrides.data())),
1087 "FFTWPlan::execute(): strides mismatch between plan and input data.");
1088 vigra_precondition((outs.stride() == TinyVectorView<int, N>(outstrides.data())),
1089 "FFTWPlan::execute(): strides mismatch between plan and output data.");
1091 detail::fftwPlanExecute(plan, ins.data(), outs.data());
1093 typedef typename MO::value_type V;
1094 if(sign == FFTW_BACKWARD)
1095 outs *= V(1.0) / Real(outs.size());
1137 template <
unsigned int N,
class Real>
1145 RArray realArray, realKernel;
1146 CArray fourierArray, fourierKernel;
1147 bool useFourierKernel;
1158 : useFourierKernel(false)
1172 template <
class C1,
class C2,
class C3>
1176 unsigned int planner_flags = FFTW_ESTIMATE)
1177 : useFourierKernel(false)
1179 init(in, kernel, out, planner_flags);
1193 template <
class C1,
class C2,
class C3>
1197 unsigned int planner_flags = FFTW_ESTIMATE)
1198 : useFourierKernel(true)
1200 init(in, kernel, out, planner_flags);
1215 template <
class C1,
class C2,
class C3>
1219 bool fourierDomainKernel,
1220 unsigned int planner_flags = FFTW_ESTIMATE)
1222 init(in, kernel, out, fourierDomainKernel, planner_flags);
1238 template <
class C1,
class C2,
class C3>
1240 bool useFourierKernel =
false,
1241 unsigned int planner_flags = FFTW_ESTIMATE)
1243 if(useFourierKernel)
1244 init(inOut, kernel, planner_flags);
1246 initFourierKernel(inOut, kernel, planner_flags);
1253 template <
class C1,
class C2,
class C3>
1257 unsigned int planner_flags = FFTW_ESTIMATE)
1259 vigra_precondition(in.
shape() == out.
shape(),
1260 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1268 template <
class C1,
class C2,
class C3>
1272 unsigned int planner_flags = FFTW_ESTIMATE)
1274 vigra_precondition(in.
shape() == out.
shape(),
1275 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1276 initFourierKernel(in.
shape(), kernel.shape(), planner_flags);
1283 template <
class C1,
class C2,
class C3>
1287 bool fourierDomainKernel,
1288 unsigned int planner_flags = FFTW_ESTIMATE)
1290 vigra_precondition(in.shape() == out.shape(),
1291 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1292 useFourierKernel = fourierDomainKernel;
1293 initComplex(in.shape(), kernel.shape(), planner_flags);
1302 template <
class C1,
class KernelIterator,
class OutIterator>
1304 KernelIterator kernels, KernelIterator kernelsEnd,
1305 OutIterator outs,
unsigned int planner_flags = FFTW_ESTIMATE)
1307 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1308 typedef typename KernelArray::value_type KernelValue;
1309 typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1310 typedef typename OutArray::value_type OutValue;
1312 bool realKernel = IsSameType<KernelValue, Real>::value;
1313 bool fourierKernel = IsSameType<KernelValue, Complex>::value;
1315 vigra_precondition(realKernel || fourierKernel,
1316 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
1317 vigra_precondition((IsSameType<OutValue, Real>::value),
1318 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
1327 initFourierKernelMany(in.
shape(),
1328 checkShapesFourier(in.
shape(), kernels, kernelsEnd, outs),
1339 template <
class C1,
class KernelIterator,
class OutIterator>
1341 KernelIterator kernels, KernelIterator kernelsEnd,
1343 bool fourierDomainKernels,
1344 unsigned int planner_flags = FFTW_ESTIMATE)
1346 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1347 typedef typename KernelArray::value_type KernelValue;
1348 typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1349 typedef typename OutArray::value_type OutValue;
1351 vigra_precondition((IsSameType<KernelValue, Complex>::value),
1352 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
1353 vigra_precondition((IsSameType<OutValue, Complex>::value),
1354 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
1356 useFourierKernel = fourierDomainKernels;
1358 Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs);
1360 CArray newFourierArray(paddedShape), newFourierKernel(paddedShape);
1362 FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags);
1363 FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags);
1365 forward_plan = fplan;
1366 backward_plan = bplan;
1367 fourierArray.
swap(newFourierArray);
1368 fourierKernel.
swap(newFourierKernel);
1371 void init(Shape inOut, Shape kernel,
1372 unsigned int planner_flags = FFTW_ESTIMATE);
1374 void initFourierKernel(Shape inOut, Shape kernel,
1375 unsigned int planner_flags = FFTW_ESTIMATE);
1377 void initComplex(Shape inOut, Shape kernel,
1378 unsigned int planner_flags = FFTW_ESTIMATE);
1380 void initMany(Shape inOut, Shape maxKernel,
1381 unsigned int planner_flags = FFTW_ESTIMATE)
1383 init(inOut, maxKernel, planner_flags);
1386 void initFourierKernelMany(Shape inOut, Shape kernels,
1387 unsigned int planner_flags = FFTW_ESTIMATE)
1389 initFourierKernel(inOut, kernels, planner_flags);
1399 template <
class C1,
class C2,
class C3>
1400 void execute(MultiArrayView<N, Real, C1> in,
1401 MultiArrayView<N, Real, C2> kernel,
1402 MultiArrayView<N, Real, C3> out);
1411 template <
class C1,
class C2,
class C3>
1412 void execute(MultiArrayView<N, Real, C1> in,
1413 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1414 MultiArrayView<N, Real, C3> out);
1423 template <
class C1,
class C2,
class C3>
1424 void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1425 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1426 MultiArrayView<N, FFTWComplex<Real>, C3> out);
1436 template <
class C1,
class KernelIterator,
class OutIterator>
1437 void executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1438 KernelIterator kernels, KernelIterator kernelsEnd,
1448 template <
class C1,
class KernelIterator,
class OutIterator>
1450 KernelIterator kernels, KernelIterator kernelsEnd,
1453 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1454 typedef typename KernelArray::value_type KernelValue;
1455 typedef typename IsSameType<KernelValue, Complex>::type UseFourierKernel;
1456 typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1457 typedef typename OutArray::value_type OutValue;
1459 bool realKernel = IsSameType<KernelValue, Real>::value;
1460 bool fourierKernel = IsSameType<KernelValue, Complex>::value;
1462 vigra_precondition(realKernel || fourierKernel,
1463 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type.");
1464 vigra_precondition((IsSameType<OutValue, Real>::value),
1465 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type.");
1467 executeManyImpl(in, kernels, kernelsEnd, outs, UseFourierKernel());
1472 template <
class KernelIterator,
class OutIterator>
1473 Shape checkShapes(Shape in,
1474 KernelIterator kernels, KernelIterator kernelsEnd,
1477 template <
class KernelIterator,
class OutIterator>
1478 Shape checkShapesFourier(Shape in,
1479 KernelIterator kernels, KernelIterator kernelsEnd,
1482 template <
class KernelIterator,
class OutIterator>
1483 Shape checkShapesComplex(Shape in,
1484 KernelIterator kernels, KernelIterator kernelsEnd,
1487 template <
class C1,
class KernelIterator,
class OutIterator>
1490 KernelIterator kernels, KernelIterator kernelsEnd,
1491 OutIterator outs, VigraFalseType );
1493 template <
class C1,
class KernelIterator,
class OutIterator>
1496 KernelIterator kernels, KernelIterator kernelsEnd,
1497 OutIterator outs, VigraTrueType );
1501 template <
unsigned int N,
class Real>
1504 unsigned int planner_flags)
1506 Shape paddedShape = fftwBestPaddedShapeR2C(in + kernel - Shape(1)),
1509 CArray newFourierArray(complexShape), newFourierKernel(complexShape);
1511 Shape realStrides = 2*newFourierArray.stride();
1513 RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data());
1514 RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data());
1516 FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags);
1517 FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags);
1519 forward_plan = fplan;
1520 backward_plan = bplan;
1521 realArray = newRealArray;
1522 realKernel = newRealKernel;
1523 fourierArray.swap(newFourierArray);
1524 fourierKernel.swap(newFourierKernel);
1525 useFourierKernel =
false;
1528 template <
unsigned int N,
class Real>
1530 FFTWConvolvePlan<N, Real>::initFourierKernel(Shape in, Shape kernel,
1531 unsigned int planner_flags)
1533 Shape complexShape = kernel,
1534 paddedShape = fftwCorrespondingShapeC2R(complexShape);
1536 for(
unsigned int k=0; k<N; ++k)
1537 vigra_precondition(in[k] <= paddedShape[k],
1538 "FFTWConvolvePlan::init(): kernel too small for given input.");
1540 CArray newFourierArray(complexShape), newFourierKernel(complexShape);
1542 Shape realStrides = 2*newFourierArray.stride();
1544 RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data());
1545 RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data());
1547 FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags);
1548 FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags);
1550 forward_plan = fplan;
1551 backward_plan = bplan;
1552 realArray = newRealArray;
1553 realKernel = newRealKernel;
1554 fourierArray.swap(newFourierArray);
1555 fourierKernel.swap(newFourierKernel);
1556 useFourierKernel =
true;
1559 template <
unsigned int N,
class Real>
1561 FFTWConvolvePlan<N, Real>::initComplex(Shape in, Shape kernel,
1562 unsigned int planner_flags)
1566 if(useFourierKernel)
1568 for(
unsigned int k=0; k<N; ++k)
1569 vigra_precondition(in[k] <= kernel[k],
1570 "FFTWConvolvePlan::init(): kernel too small for given input.");
1572 paddedShape = kernel;
1576 paddedShape = fftwBestPaddedShape(in + kernel - Shape(1));
1579 CArray newFourierArray(paddedShape), newFourierKernel(paddedShape);
1581 FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags);
1582 FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags);
1584 forward_plan = fplan;
1585 backward_plan = bplan;
1586 fourierArray.swap(newFourierArray);
1587 fourierKernel.swap(newFourierKernel);
1590 #ifndef DOXYGEN // doxygen documents these functions as free functions
1592 template <
unsigned int N,
class Real>
1593 template <
class C1,
class C2,
class C3>
1596 MultiArrayView<N, Real, C2> kernel,
1597 MultiArrayView<N, Real, C3> out)
1599 vigra_precondition(!useFourierKernel,
1600 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel.");
1602 vigra_precondition(in.shape() == out.shape(),
1603 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1605 Shape paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernel.shape() - Shape(1)),
1606 diff = paddedShape - in.shape(),
1608 right = in.shape() + left;
1610 vigra_precondition(paddedShape == realArray.shape(),
1611 "FFTWConvolvePlan::execute(): shape mismatch between input and plan.");
1613 detail::fftEmbedArray(in, realArray);
1614 forward_plan.execute(realArray, fourierArray);
1616 detail::fftEmbedKernel(kernel, realKernel);
1617 forward_plan.execute(realKernel, fourierKernel);
1619 fourierArray *= fourierKernel;
1621 backward_plan.execute(fourierArray, realArray);
1623 out = realArray.subarray(left, right);
1626 template <
unsigned int N,
class Real>
1627 template <
class C1,
class C2,
class C3>
1630 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1631 MultiArrayView<N, Real, C3> out)
1633 vigra_precondition(useFourierKernel,
1634 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel.");
1636 vigra_precondition(in.shape() == out.shape(),
1637 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1639 vigra_precondition(kernel.shape() == fourierArray.shape(),
1640 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan.");
1642 Shape paddedShape = fftwCorrespondingShapeC2R(kernel.shape(),
odd(in.shape(0))),
1643 diff = paddedShape - in.shape(),
1645 right = in.shape() + left;
1647 vigra_precondition(paddedShape == realArray.shape(),
1648 "FFTWConvolvePlan::execute(): shape mismatch between input and plan.");
1650 detail::fftEmbedArray(in, realArray);
1651 forward_plan.execute(realArray, fourierArray);
1653 fourierKernel = kernel;
1654 moveDCToHalfspaceUpperLeft(fourierKernel);
1656 fourierArray *= fourierKernel;
1658 backward_plan.execute(fourierArray, realArray);
1660 out = realArray.subarray(left, right);
1663 template <
unsigned int N,
class Real>
1664 template <
class C1,
class C2,
class C3>
1667 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1668 MultiArrayView<N, FFTWComplex<Real>, C3> out)
1670 vigra_precondition(in.shape() == out.shape(),
1671 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1673 Shape paddedShape = fourierArray.shape(),
1674 diff = paddedShape - in.shape(),
1676 right = in.shape() + left;
1678 if(useFourierKernel)
1680 vigra_precondition(kernel.shape() == fourierArray.shape(),
1681 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan.");
1683 fourierKernel = kernel;
1688 detail::fftEmbedKernel(kernel, fourierKernel);
1689 forward_plan.execute(fourierKernel, fourierKernel);
1692 detail::fftEmbedArray(in, fourierArray);
1693 forward_plan.execute(fourierArray, fourierArray);
1695 fourierArray *= fourierKernel;
1697 backward_plan.execute(fourierArray, fourierArray);
1699 out = fourierArray.subarray(left, right);
1702 template <
unsigned int N,
class Real>
1703 template <
class C1,
class KernelIterator,
class OutIterator>
1705 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
1706 KernelIterator kernels, KernelIterator kernelsEnd,
1707 OutIterator outs, VigraFalseType )
1709 vigra_precondition(!useFourierKernel,
1710 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel.");
1712 Shape kernelMax = checkShapes(in.shape(), kernels, kernelsEnd, outs),
1713 paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernelMax - Shape(1)),
1714 diff = paddedShape - in.shape(),
1716 right = in.shape() + left;
1718 vigra_precondition(paddedShape == realArray.shape(),
1719 "FFTWConvolvePlan::executeMany(): shape mismatch between input and plan.");
1721 detail::fftEmbedArray(in, realArray);
1722 forward_plan.execute(realArray, fourierArray);
1724 for(; kernels != kernelsEnd; ++kernels, ++outs)
1726 detail::fftEmbedKernel(*kernels, realKernel);
1727 forward_plan.execute(realKernel, fourierKernel);
1729 fourierKernel *= fourierArray;
1731 backward_plan.execute(fourierKernel, realKernel);
1733 *outs = realKernel.subarray(left, right);
1737 template <
unsigned int N,
class Real>
1738 template <
class C1,
class KernelIterator,
class OutIterator>
1740 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
1741 KernelIterator kernels, KernelIterator kernelsEnd,
1742 OutIterator outs, VigraTrueType )
1744 vigra_precondition(useFourierKernel,
1745 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel.");
1747 Shape complexShape = checkShapesFourier(in.shape(), kernels, kernelsEnd, outs),
1748 paddedShape = fftwCorrespondingShapeC2R(complexShape,
odd(in.shape(0))),
1749 diff = paddedShape - in.shape(),
1751 right = in.shape() + left;
1753 vigra_precondition(complexShape == fourierArray.shape(),
1754 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between kernels and plan.");
1756 vigra_precondition(paddedShape == realArray.shape(),
1757 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between input and plan.");
1759 detail::fftEmbedArray(in, realArray);
1760 forward_plan.execute(realArray, fourierArray);
1762 for(; kernels != kernelsEnd; ++kernels, ++outs)
1764 fourierKernel = *kernels;
1765 moveDCToHalfspaceUpperLeft(fourierKernel);
1766 fourierKernel *= fourierArray;
1768 backward_plan.execute(fourierKernel, realKernel);
1770 *outs = realKernel.subarray(left, right);
1774 template <
unsigned int N,
class Real>
1775 template <
class C1,
class KernelIterator,
class OutIterator>
1778 KernelIterator kernels, KernelIterator kernelsEnd,
1781 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1782 typedef typename KernelArray::value_type KernelValue;
1783 typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1784 typedef typename OutArray::value_type OutValue;
1786 vigra_precondition((IsSameType<KernelValue, Complex>::value),
1787 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type.");
1788 vigra_precondition((IsSameType<OutValue, Complex>::value),
1789 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type.");
1791 Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs),
1792 diff = paddedShape - in.shape(),
1794 right = in.shape() + left;
1796 detail::fftEmbedArray(in, fourierArray);
1797 forward_plan.execute(fourierArray, fourierArray);
1799 for(; kernels != kernelsEnd; ++kernels, ++outs)
1801 if(useFourierKernel)
1803 fourierKernel = *kernels;
1808 detail::fftEmbedKernel(*kernels, fourierKernel);
1809 forward_plan.execute(fourierKernel, fourierKernel);
1812 fourierKernel *= fourierArray;
1814 backward_plan.execute(fourierKernel, fourierKernel);
1816 *outs = fourierKernel.subarray(left, right);
1822 template <
unsigned int N,
class Real>
1823 template <
class KernelIterator,
class OutIterator>
1824 typename FFTWConvolvePlan<N, Real>::Shape
1825 FFTWConvolvePlan<N, Real>::checkShapes(Shape in,
1826 KernelIterator kernels, KernelIterator kernelsEnd,
1829 vigra_precondition(kernels != kernelsEnd,
1830 "FFTWConvolvePlan::checkShapes(): empty kernel sequence.");
1833 for(; kernels != kernelsEnd; ++kernels, ++outs)
1835 vigra_precondition(in == outs->shape(),
1836 "FFTWConvolvePlan::checkShapes(): shape mismatch between input and (one) output.");
1837 kernelMax =
max(kernelMax, kernels->shape());
1839 vigra_precondition(
prod(kernelMax) > 0,
1840 "FFTWConvolvePlan::checkShapes(): all kernels have size 0.");
1844 template <
unsigned int N,
class Real>
1845 template <
class KernelIterator,
class OutIterator>
1846 typename FFTWConvolvePlan<N, Real>::Shape
1847 FFTWConvolvePlan<N, Real>::checkShapesFourier(Shape in,
1848 KernelIterator kernels, KernelIterator kernelsEnd,
1851 vigra_precondition(kernels != kernelsEnd,
1852 "FFTWConvolvePlan::checkShapesFourier(): empty kernel sequence.");
1854 Shape complexShape = kernels->shape(),
1855 paddedShape = fftwCorrespondingShapeC2R(complexShape);
1857 for(
unsigned int k=0; k<N; ++k)
1858 vigra_precondition(in[k] <= paddedShape[k],
1859 "FFTWConvolvePlan::checkShapesFourier(): kernels too small for given input.");
1861 for(; kernels != kernelsEnd; ++kernels, ++outs)
1863 vigra_precondition(in == outs->shape(),
1864 "FFTWConvolvePlan::checkShapesFourier(): shape mismatch between input and (one) output.");
1865 vigra_precondition(complexShape == kernels->shape(),
1866 "FFTWConvolvePlan::checkShapesFourier(): all kernels must have the same size.");
1868 return complexShape;
1871 template <
unsigned int N,
class Real>
1872 template <
class KernelIterator,
class OutIterator>
1873 typename FFTWConvolvePlan<N, Real>::Shape
1874 FFTWConvolvePlan<N, Real>::checkShapesComplex(Shape in,
1875 KernelIterator kernels, KernelIterator kernelsEnd,
1878 vigra_precondition(kernels != kernelsEnd,
1879 "FFTWConvolvePlan::checkShapesComplex(): empty kernel sequence.");
1881 Shape kernelShape = kernels->shape();
1882 for(; kernels != kernelsEnd; ++kernels, ++outs)
1884 vigra_precondition(in == outs->shape(),
1885 "FFTWConvolvePlan::checkShapesComplex(): shape mismatch between input and (one) output.");
1886 if(useFourierKernel)
1888 vigra_precondition(kernelShape == kernels->shape(),
1889 "FFTWConvolvePlan::checkShapesComplex(): Fourier domain kernels must have identical size.");
1893 kernelShape =
max(kernelShape, kernels->shape());
1896 vigra_precondition(
prod(kernelShape) > 0,
1897 "FFTWConvolvePlan::checkShapesComplex(): all kernels have size 0.");
1899 if(useFourierKernel)
1901 for(
unsigned int k=0; k<N; ++k)
1902 vigra_precondition(in[k] <= kernelShape[k],
1903 "FFTWConvolvePlan::checkShapesComplex(): kernels too small for given input.");
1908 return fftwBestPaddedShape(in + kernelShape - Shape(1));
1919 template <
unsigned int N,
class Real,
class C1,
class C2>
1922 MultiArrayView<N, FFTWComplex<Real>, C2> out)
1924 FFTWPlan<N, Real>(in, out, FFTW_FORWARD).execute(in, out);
1927 template <
unsigned int N,
class Real,
class C1,
class C2>
1930 MultiArrayView<N, FFTWComplex<Real>, C2> out)
1932 FFTWPlan<N, Real>(in, out, FFTW_BACKWARD).execute(in, out);
1935 template <
unsigned int N,
class Real,
class C1,
class C2>
1938 MultiArrayView<N, FFTWComplex<Real>, C2> out)
1940 if(in.shape() == out.shape())
1944 FFTWPlan<N, Real>(out, out, FFTW_FORWARD).execute(out, out);
1948 FFTWPlan<N, Real>(in, out).execute(in, out);
1951 vigra_precondition(
false,
1952 "fourierTransform(): shape mismatch between input and output.");
1955 template <
unsigned int N,
class Real,
class C1,
class C2>
1958 MultiArrayView<N, Real, C2> out)
1961 "fourierTransformInverse(): shape mismatch between input and output.");
1962 FFTWPlan<N, Real>(in, out).execute(in, out);
2154 doxygen_overloaded_function(template <...>
void convolveFFT)
2156 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2159 MultiArrayView<N, Real, C2> kernel,
2160 MultiArrayView<N, Real, C3> out)
2162 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out);
2165 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2168 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2169 MultiArrayView<N, Real, C3> out)
2171 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out);
2180 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2183 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2184 MultiArrayView<N, FFTWComplex<Real>, C3> out,
2185 bool fourierDomainKernel)
2187 FFTWConvolvePlan<N, Real>(in, kernel, out, fourierDomainKernel).execute(in, kernel, out);
2196 template <
unsigned int N,
class Real,
class C1,
2197 class KernelIterator,
class OutIterator>
2200 KernelIterator kernels, KernelIterator kernelsEnd,
2203 FFTWConvolvePlan<N, Real> plan;
2204 plan.initMany(in, kernels, kernelsEnd, outs);
2205 plan.executeMany(in, kernels, kernelsEnd, outs);
2214 template <
unsigned int N,
class Real,
class C1,
2215 class KernelIterator,
class OutIterator>
2218 KernelIterator kernels, KernelIterator kernelsEnd,
2220 bool fourierDomainKernel)
2222 FFTWConvolvePlan<N, Real> plan;
2223 plan.initMany(in, kernels, kernelsEnd, outs, fourierDomainKernel);
2224 plan.executeMany(in, kernels, kernelsEnd, outs);
2231 #endif // VIGRA_MULTI_FFT_HXX