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 struct FFTWPaddingSize
370 static const int size = 1330;
371 static const int evenSize = 1063;
372 static int goodSizes[size];
373 static int goodEvenSizes[evenSize];
375 static inline int find(
int s)
377 if(s <= 0 || s >= goodSizes[size-1])
380 int * upperBound = std::upper_bound(goodSizes, goodSizes+size, s);
381 if(upperBound > goodSizes && upperBound[-1] == s)
387 static inline int findEven(
int s)
389 if(s <= 0 || s >= goodEvenSizes[evenSize-1])
392 int * upperBound = std::upper_bound(goodEvenSizes, goodEvenSizes+evenSize, s);
393 if(upperBound > goodEvenSizes && upperBound[-1] == s)
405 int FFTWPaddingSize<DUMMY>::goodSizes[size] = {
406 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
407 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 35, 36, 39, 40, 42, 44, 45, 48,
408 49, 50, 52, 54, 55, 56, 60, 63, 64, 65, 66, 70, 72, 75, 77, 78, 80, 81,
409 84, 88, 90, 91, 96, 98, 99, 100, 104, 105, 108, 110, 112, 117, 120, 125,
410 126, 128, 130, 132, 135, 140, 144, 147, 150, 154, 156, 160, 162, 165,
411 168, 175, 176, 180, 182, 189, 192, 195, 196, 198, 200, 208, 210, 216,
412 220, 224, 225, 231, 234, 240, 243, 245, 250, 252, 256, 260, 264, 270,
413 273, 275, 280, 288, 294, 297, 300, 308, 312, 315, 320, 324, 325, 330,
414 336, 343, 350, 351, 352, 360, 364, 375, 378, 384, 385, 390, 392, 396,
415 400, 405, 416, 420, 432, 440, 441, 448, 450, 455, 462, 468, 480, 486,
416 490, 495, 500, 504, 512, 520, 525, 528, 539, 540, 546, 550, 560, 567,
417 576, 585, 588, 594, 600, 616, 624, 625, 630, 637, 640, 648, 650, 660,
418 672, 675, 686, 693, 700, 702, 704, 720, 728, 729, 735, 750, 756, 768,
419 770, 780, 784, 792, 800, 810, 819, 825, 832, 840, 864, 875, 880, 882,
420 891, 896, 900, 910, 924, 936, 945, 960, 972, 975, 980, 990, 1000, 1008,
421 1024, 1029, 1040, 1050, 1053, 1056, 1078, 1080, 1092, 1100, 1120, 1125,
422 1134, 1152, 1155, 1170, 1176, 1188, 1200, 1215, 1225, 1232, 1248, 1250,
423 1260, 1274, 1280, 1296, 1300, 1320, 1323, 1344, 1350, 1365, 1372, 1375,
424 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470, 1485, 1500, 1512, 1536,
425 1540, 1560, 1568, 1575, 1584, 1600, 1617, 1620, 1625, 1638, 1650, 1664,
426 1680, 1701, 1715, 1728, 1750, 1755, 1760, 1764, 1782, 1792, 1800, 1820,
427 1848, 1872, 1875, 1890, 1911, 1920, 1925, 1944, 1950, 1960, 1980, 2000,
428 2016, 2025, 2048, 2058, 2079, 2080, 2100, 2106, 2112, 2156, 2160, 2184,
429 2187, 2200, 2205, 2240, 2250, 2268, 2275, 2304, 2310, 2340, 2352, 2376,
430 2400, 2401, 2430, 2450, 2457, 2464, 2475, 2496, 2500, 2520, 2548, 2560,
431 2592, 2600, 2625, 2640, 2646, 2673, 2688, 2695, 2700, 2730, 2744, 2750,
432 2772, 2800, 2808, 2816, 2835, 2880, 2912, 2916, 2925, 2940, 2970, 3000,
433 3024, 3072, 3080, 3087, 3120, 3125, 3136, 3150, 3159, 3168, 3185, 3200,
434 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3375, 3402, 3430, 3456, 3465,
435 3500, 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3645, 3675, 3696, 3744,
436 3750, 3773, 3780, 3822, 3840, 3850, 3888, 3900, 3920, 3960, 3969, 4000,
437 4032, 4050, 4095, 4096, 4116, 4125, 4158, 4160, 4200, 4212, 4224, 4312,
438 4320, 4368, 4374, 4375, 4400, 4410, 4455, 4459, 4480, 4500, 4536, 4550,
439 4608, 4620, 4680, 4704, 4725, 4752, 4800, 4802, 4851, 4860, 4875, 4900,
440 4914, 4928, 4950, 4992, 5000, 5040, 5096, 5103, 5120, 5145, 5184, 5200,
441 5250, 5265, 5280, 5292, 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544,
442 5600, 5616, 5625, 5632, 5670, 5733, 5760, 5775, 5824, 5832, 5850, 5880,
443 5940, 6000, 6048, 6075, 6125, 6144, 6160, 6174, 6237, 6240, 6250, 6272,
444 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6561, 6600, 6615,
445 6656, 6720, 6750, 6804, 6825, 6860, 6875, 6912, 6930, 7000, 7020, 7040,
446 7056, 7128, 7168, 7200, 7203, 7280, 7290, 7350, 7371, 7392, 7425, 7488,
447 7500, 7546, 7560, 7644, 7680, 7700, 7776, 7800, 7840, 7875, 7920, 7938,
448 8000, 8019, 8064, 8085, 8100, 8125, 8190, 8192, 8232, 8250, 8316, 8320,
449 8400, 8424, 8448, 8505, 8575, 8624, 8640, 8736, 8748, 8750, 8775, 8800,
450 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9261, 9360, 9375,
451 9408, 9450, 9477, 9504, 9555, 9600, 9604, 9625, 9702, 9720, 9750, 9800,
452 9828, 9856, 9900, 9984, 10000, 10080, 10125, 10192, 10206, 10240, 10290,
453 10368, 10395, 10400, 10500, 10530, 10560, 10584, 10692, 10752, 10780,
454 10800, 10920, 10935, 10976, 11000, 11025, 11088, 11200, 11232, 11250,
455 11264, 11319, 11340, 11375, 11466, 11520, 11550, 11648, 11664, 11700,
456 11760, 11880, 11907, 12000, 12005, 12096, 12150, 12250, 12285, 12288,
457 12320, 12348, 12375, 12474, 12480, 12500, 12544, 12600, 12636, 12672,
458 12740, 12800, 12936, 12960, 13000, 13104, 13122, 13125, 13200, 13230,
459 13312, 13365, 13377, 13440, 13475, 13500, 13608, 13650, 13720, 13750,
460 13824, 13860, 14000, 14040, 14080, 14112, 14175, 14256, 14336, 14400,
461 14406, 14553, 14560, 14580, 14625, 14700, 14742, 14784, 14850, 14976,
462 15000, 15092, 15120, 15288, 15309, 15360, 15400, 15435, 15552, 15600,
463 15625, 15680, 15750, 15795, 15840, 15876, 15925, 16000, 16038, 16128,
464 16170, 16200, 16250, 16380, 16384, 16464, 16500, 16632, 16640, 16800,
465 16807, 16848, 16875, 16896, 17010, 17150, 17199, 17248, 17280, 17325,
466 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000,
467 18144, 18200, 18225, 18375, 18432, 18480, 18522, 18711, 18720, 18750,
468 18816, 18865, 18900, 18954, 19008, 19110, 19200, 19208, 19250, 19404,
469 19440, 19500, 19600, 19656, 19683, 19712, 19800, 19845, 19968, 20000,
470 20160, 20250, 20384, 20412, 20475, 20480, 20580, 20625, 20736, 20790,
471 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560, 21600, 21609,
472 21840, 21870, 21875, 21952, 22000, 22050, 22113, 22176, 22275, 22295,
473 22400, 22464, 22500, 22528, 22638, 22680, 22750, 22932, 23040, 23100,
474 23296, 23328, 23400, 23520, 23625, 23760, 23814, 24000, 24010, 24057,
475 24192, 24255, 24300, 24375, 24500, 24570, 24576, 24640, 24696, 24750,
476 24948, 24960, 25000, 25088, 25200, 25272, 25344, 25480, 25515, 25600,
477 25725, 25872, 25920, 26000, 26208, 26244, 26250, 26325, 26400, 26411,
478 26460, 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440,
479 27500, 27648, 27720, 27783, 28000, 28080, 28125, 28160, 28224, 28350,
480 28431, 28512, 28665, 28672, 28800, 28812, 28875, 29106, 29120, 29160,
481 29250, 29400, 29484, 29568, 29700, 29952, 30000, 30184, 30240, 30375,
482 30576, 30618, 30625, 30720, 30800, 30870, 31104, 31185, 31200, 31213,
483 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256,
484 32340, 32400, 32500, 32760, 32768, 32805, 32928, 33000, 33075, 33264,
485 33280, 33600, 33614, 33696, 33750, 33792, 33957, 34020, 34125, 34300,
486 34375, 34398, 34496, 34560, 34650, 34944, 34992, 35000, 35100, 35200,
487 35280, 35640, 35672, 35721, 35840, 36000, 36015, 36288, 36400, 36450,
488 36750, 36855, 36864, 36960, 37044, 37125, 37422, 37440, 37500, 37632,
489 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808, 38880,
490 39000, 39200, 39312, 39366, 39375, 39424, 39600, 39690, 39936, 40000,
491 40095, 40131, 40320, 40425, 40500, 40625, 40768, 40824, 40950, 40960,
492 41160, 41250, 41472, 41580, 41600, 42000, 42120, 42240, 42336, 42525,
493 42768, 42875, 43008, 43120, 43200, 43218, 43659, 43680, 43740, 43750,
494 43875, 43904, 44000, 44100, 44226, 44352, 44550, 44590, 44800, 44928,
495 45000, 45056, 45276, 45360, 45500, 45864, 45927, 46080, 46200, 46305,
496 46592, 46656, 46800, 46875, 47040, 47250, 47385, 47520, 47628, 47775,
497 48000, 48020, 48114, 48125, 48384, 48510, 48600, 48750, 49000, 49140,
498 49152, 49280, 49392, 49500, 49896, 49920, 50000, 50176, 50400, 50421,
499 50544, 50625, 50688, 50960, 51030, 51200, 51450, 51597, 51744, 51840,
500 51975, 52000, 52416, 52488, 52500, 52650, 52800, 52822, 52920, 53248,
501 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54675, 54880, 55000,
502 55125, 55296, 55440, 55566, 56000, 56133, 56160, 56250, 56320, 56448,
503 56595, 56700, 56862, 56875, 57024, 57330, 57344, 57600, 57624, 57750,
504 58212, 58240, 58320, 58500, 58800, 58968, 59049, 59136, 59400, 59535,
505 59904, 60000, 60025, 60368, 60480, 60750, 61152, 61236, 61250, 61425,
506 61440, 61600, 61740, 61875, 62208, 62370, 62400, 62426, 62500, 62720,
507 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512, 64680, 64800,
508 64827, 65000, 65520, 65536, 65610, 65625, 65856, 66000, 66150, 66339,
509 66528, 66560, 66825, 66885, 67200, 67228, 67375, 67392, 67500, 67584,
510 67914, 68040, 68250, 68600, 68750, 68796, 68992, 69120, 69300, 69888,
511 69984, 70000, 70200, 70400, 70560, 70875, 71280, 71344, 71442, 71680,
512 72000, 72030, 72171, 72576, 72765, 72800, 72900, 73125, 73500, 73710,
513 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264, 75460, 75600,
514 75816, 76032, 76440, 76545, 76800, 76832, 77000, 77175, 77616, 77760,
515 78000, 78125, 78400, 78624, 78732, 78750, 78848, 78975, 79200, 79233,
516 79380, 79625, 79872, 80000, 80190, 80262, 80640, 80850, 81000, 81250,
517 81536, 81648, 81900, 81920, 82320, 82500, 82944, 83160, 83200, 83349,
518 84000, 84035, 84240, 84375, 84480, 84672, 85050, 85293, 85536, 85750,
519 85995, 86016, 86240, 86400, 86436, 86625, 87318, 87360, 87480, 87500,
520 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856,
521 90000, 90112, 90552, 90720, 91000, 91125, 91728, 91854, 91875, 92160,
522 92400, 92610, 93184, 93312, 93555, 93600, 93639, 93750, 94080, 94325,
523 94500, 94770, 95040, 95256, 95550, 96000, 96040, 96228, 96250, 96768,
524 97020, 97200, 97500, 98000, 98280, 98304, 98415, 98560, 98784, 99000,
529 int FFTWPaddingSize<DUMMY>::goodEvenSizes[evenSize] = {
530 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22,
531 24, 26, 28, 30, 32, 36, 40, 42, 44, 48, 50, 52, 54, 56, 60, 64, 66, 70,
532 72, 78, 80, 84, 88, 90, 96, 98, 100, 104, 108, 110, 112, 120, 126, 128,
533 130, 132, 140, 144, 150, 154, 156, 160, 162, 168, 176, 180, 182, 192,
534 196, 198, 200, 208, 210, 216, 220, 224, 234, 240, 250, 252, 256, 260,
535 264, 270, 280, 288, 294, 300, 308, 312, 320, 324, 330, 336, 350, 352,
536 360, 364, 378, 384, 390, 392, 396, 400, 416, 420, 432, 440, 448, 450,
537 462, 468, 480, 486, 490, 500, 504, 512, 520, 528, 540, 546, 550, 560,
538 576, 588, 594, 600, 616, 624, 630, 640, 648, 650, 660, 672, 686, 700,
539 702, 704, 720, 728, 750, 756, 768, 770, 780, 784, 792, 800, 810, 832,
540 840, 864, 880, 882, 896, 900, 910, 924, 936, 960, 972, 980, 990, 1000,
541 1008, 1024, 1040, 1050, 1056, 1078, 1080, 1092, 1100, 1120, 1134, 1152,
542 1170, 1176, 1188, 1200, 1232, 1248, 1250, 1260, 1274, 1280, 1296, 1300,
543 1320, 1344, 1350, 1372, 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470,
544 1500, 1512, 1536, 1540, 1560, 1568, 1584, 1600, 1620, 1638, 1650, 1664,
545 1680, 1728, 1750, 1760, 1764, 1782, 1792, 1800, 1820, 1848, 1872, 1890,
546 1920, 1944, 1950, 1960, 1980, 2000, 2016, 2048, 2058, 2080, 2100, 2106,
547 2112, 2156, 2160, 2184, 2200, 2240, 2250, 2268, 2304, 2310, 2340, 2352,
548 2376, 2400, 2430, 2450, 2464, 2496, 2500, 2520, 2548, 2560, 2592, 2600,
549 2640, 2646, 2688, 2700, 2730, 2744, 2750, 2772, 2800, 2808, 2816, 2880,
550 2912, 2916, 2940, 2970, 3000, 3024, 3072, 3080, 3120, 3136, 3150, 3168,
551 3200, 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3402, 3430, 3456, 3500,
552 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3696, 3744, 3750, 3780, 3822,
553 3840, 3850, 3888, 3900, 3920, 3960, 4000, 4032, 4050, 4096, 4116, 4158,
554 4160, 4200, 4212, 4224, 4312, 4320, 4368, 4374, 4400, 4410, 4480, 4500,
555 4536, 4550, 4608, 4620, 4680, 4704, 4752, 4800, 4802, 4860, 4900, 4914,
556 4928, 4950, 4992, 5000, 5040, 5096, 5120, 5184, 5200, 5250, 5280, 5292,
557 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544, 5600, 5616, 5632, 5670,
558 5760, 5824, 5832, 5850, 5880, 5940, 6000, 6048, 6144, 6160, 6174, 6240,
559 6250, 6272, 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6600,
560 6656, 6720, 6750, 6804, 6860, 6912, 6930, 7000, 7020, 7040, 7056, 7128,
561 7168, 7200, 7280, 7290, 7350, 7392, 7488, 7500, 7546, 7560, 7644, 7680,
562 7700, 7776, 7800, 7840, 7920, 7938, 8000, 8064, 8100, 8190, 8192, 8232,
563 8250, 8316, 8320, 8400, 8424, 8448, 8624, 8640, 8736, 8748, 8750, 8800,
564 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9360, 9408, 9450,
565 9504, 9600, 9604, 9702, 9720, 9750, 9800, 9828, 9856, 9900, 9984, 10000,
566 10080, 10192, 10206, 10240, 10290, 10368, 10400, 10500, 10530, 10560,
567 10584, 10692, 10752, 10780, 10800, 10920, 10976, 11000, 11088, 11200,
568 11232, 11250, 11264, 11340, 11466, 11520, 11550, 11648, 11664, 11700,
569 11760, 11880, 12000, 12096, 12150, 12250, 12288, 12320, 12348, 12474,
570 12480, 12500, 12544, 12600, 12636, 12672, 12740, 12800, 12936, 12960,
571 13000, 13104, 13122, 13200, 13230, 13312, 13440, 13500, 13608, 13650,
572 13720, 13750, 13824, 13860, 14000, 14040, 14080, 14112, 14256, 14336,
573 14400, 14406, 14560, 14580, 14700, 14742, 14784, 14850, 14976, 15000,
574 15092, 15120, 15288, 15360, 15400, 15552, 15600, 15680, 15750, 15840,
575 15876, 16000, 16038, 16128, 16170, 16200, 16250, 16380, 16384, 16464,
576 16500, 16632, 16640, 16800, 16848, 16896, 17010, 17150, 17248, 17280,
577 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000,
578 18144, 18200, 18432, 18480, 18522, 18720, 18750, 18816, 18900, 18954,
579 19008, 19110, 19200, 19208, 19250, 19404, 19440, 19500, 19600, 19656,
580 19712, 19800, 19968, 20000, 20160, 20250, 20384, 20412, 20480, 20580,
581 20736, 20790, 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560,
582 21600, 21840, 21870, 21952, 22000, 22050, 22176, 22400, 22464, 22500,
583 22528, 22638, 22680, 22750, 22932, 23040, 23100, 23296, 23328, 23400,
584 23520, 23760, 23814, 24000, 24010, 24192, 24300, 24500, 24570, 24576,
585 24640, 24696, 24750, 24948, 24960, 25000, 25088, 25200, 25272, 25344,
586 25480, 25600, 25872, 25920, 26000, 26208, 26244, 26250, 26400, 26460,
587 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440, 27500,
588 27648, 27720, 28000, 28080, 28160, 28224, 28350, 28512, 28672, 28800,
589 28812, 29106, 29120, 29160, 29250, 29400, 29484, 29568, 29700, 29952,
590 30000, 30184, 30240, 30576, 30618, 30720, 30800, 30870, 31104, 31200,
591 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256,
592 32340, 32400, 32500, 32760, 32768, 32928, 33000, 33264, 33280, 33600,
593 33614, 33696, 33750, 33792, 34020, 34300, 34398, 34496, 34560, 34650,
594 34944, 34992, 35000, 35100, 35200, 35280, 35640, 35672, 35840, 36000,
595 36288, 36400, 36450, 36750, 36864, 36960, 37044, 37422, 37440, 37500,
596 37632, 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808,
597 38880, 39000, 39200, 39312, 39366, 39424, 39600, 39690, 39936, 40000,
598 40320, 40500, 40768, 40824, 40950, 40960, 41160, 41250, 41472, 41580,
599 41600, 42000, 42120, 42240, 42336, 42768, 43008, 43120, 43200, 43218,
600 43680, 43740, 43750, 43904, 44000, 44100, 44226, 44352, 44550, 44590,
601 44800, 44928, 45000, 45056, 45276, 45360, 45500, 45864, 46080, 46200,
602 46592, 46656, 46800, 47040, 47250, 47520, 47628, 48000, 48020, 48114,
603 48384, 48510, 48600, 48750, 49000, 49140, 49152, 49280, 49392, 49500,
604 49896, 49920, 50000, 50176, 50400, 50544, 50688, 50960, 51030, 51200,
605 51450, 51744, 51840, 52000, 52416, 52488, 52500, 52650, 52800, 52822,
606 52920, 53248, 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54880,
607 55000, 55296, 55440, 55566, 56000, 56160, 56250, 56320, 56448, 56700,
608 56862, 57024, 57330, 57344, 57600, 57624, 57750, 58212, 58240, 58320,
609 58500, 58800, 58968, 59136, 59400, 59904, 60000, 60368, 60480, 60750,
610 61152, 61236, 61250, 61440, 61600, 61740, 62208, 62370, 62400, 62426,
611 62500, 62720, 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512,
612 64680, 64800, 65000, 65520, 65536, 65610, 65856, 66000, 66150, 66528,
613 66560, 67200, 67228, 67392, 67500, 67584, 67914, 68040, 68250, 68600,
614 68750, 68796, 68992, 69120, 69300, 69888, 69984, 70000, 70200, 70400,
615 70560, 71280, 71344, 71442, 71680, 72000, 72030, 72576, 72800, 72900,
616 73500, 73710, 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264,
617 75460, 75600, 75816, 76032, 76440, 76800, 76832, 77000, 77616, 77760,
618 78000, 78400, 78624, 78732, 78750, 78848, 79200, 79380, 79872, 80000,
619 80190, 80262, 80640, 80850, 81000, 81250, 81536, 81648, 81900, 81920,
620 82320, 82500, 82944, 83160, 83200, 84000, 84240, 84480, 84672, 85050,
621 85536, 85750, 86016, 86240, 86400, 86436, 87318, 87360, 87480, 87500,
622 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856,
623 90000, 90112, 90552, 90720, 91000, 91728, 91854, 92160, 92400, 92610,
624 93184, 93312, 93600, 93750, 94080, 94500, 94770, 95040, 95256, 95550,
625 96000, 96040, 96228, 96250, 96768, 97020, 97200, 97500, 98000, 98280,
626 98304, 98560, 98784, 99000, 99792, 99840
630 struct FFTEmbedKernel
632 template <
unsigned int N,
class Real,
class C,
class Shape>
634 exec(MultiArrayView<N, Real, C> & out, Shape
const & kernelShape,
635 Shape & srcPoint, Shape & destPoint,
bool copyIt)
637 for(srcPoint[M]=0; srcPoint[M]<kernelShape[M]; ++srcPoint[M])
639 if(srcPoint[M] < (kernelShape[M] + 1) / 2)
641 destPoint[M] = srcPoint[M];
645 destPoint[M] = srcPoint[M] + out.shape(M) - kernelShape[M];
648 FFTEmbedKernel<M-1>::exec(out, kernelShape, srcPoint, destPoint, copyIt);
654 struct FFTEmbedKernel<0>
656 template <
unsigned int N,
class Real,
class C,
class Shape>
658 exec(MultiArrayView<N, Real, C> & out, Shape
const & kernelShape,
659 Shape & srcPoint, Shape & destPoint,
bool copyIt)
661 for(srcPoint[0]=0; srcPoint[0]<kernelShape[0]; ++srcPoint[0])
663 if(srcPoint[0] < (kernelShape[0] + 1) / 2)
665 destPoint[0] = srcPoint[0];
669 destPoint[0] = srcPoint[0] + out.shape(0) - kernelShape[0];
674 out[destPoint] = out[srcPoint];
681 template <
unsigned int N,
class Real,
class C1,
class C2>
683 fftEmbedKernel(MultiArrayView<N, Real, C1> kernel,
684 MultiArrayView<N, Real, C2> out,
687 typedef typename MultiArrayShape<N>::type Shape;
689 MultiArrayView<N, Real, C2> kout = out.subarray(Shape(), kernel.shape());
697 Shape srcPoint, destPoint;
698 FFTEmbedKernel<(int)N-1>::exec(out, kernel.shape(), srcPoint, destPoint,
false);
701 template <
unsigned int N,
class Real,
class C1,
class C2>
703 fftEmbedArray(MultiArrayView<N, Real, C1> in,
704 MultiArrayView<N, Real, C2> out)
706 typedef typename MultiArrayShape<N>::type Shape;
708 Shape diff = out.shape() - in.shape(),
710 rightDiff = diff - leftDiff,
711 right = in.shape() + leftDiff;
713 out.subarray(leftDiff, right) = in;
715 typedef typename MultiArrayView<N, Real, C2>::traverser Traverser;
716 typedef MultiArrayNavigator<Traverser, N> Navigator;
717 typedef typename Navigator::iterator Iterator;
719 for(
unsigned int d = 0; d < N; ++d)
721 Navigator nav(out.traverser_begin(), out.shape(), d);
723 for( ; nav.hasMore(); nav++ )
725 Iterator i = nav.begin();
726 for(
int k=1; k<=leftDiff[d]; ++k)
727 i[leftDiff[d] - k] = i[leftDiff[d] + k];
728 for(
int k=0; k<rightDiff[d]; ++k)
729 i[right[d] + k] = i[right[d] - k - 2];
736 template <
class T,
int N>
738 fftwBestPaddedShape(TinyVector<T, N> shape)
740 for(
unsigned int k=0; k<N; ++k)
741 shape[k] = detail::FFTWPaddingSize<0>::find(shape[k]);
745 template <
class T,
int N>
747 fftwBestPaddedShapeR2C(TinyVector<T, N> shape)
749 shape[0] = detail::FFTWPaddingSize<0>::findEven(shape[0]);
750 for(
unsigned int k=1; k<N; ++k)
751 shape[k] = detail::FFTWPaddingSize<0>::find(shape[k]);
766 template <
class T,
int N>
770 shape[0] = shape[0] / 2 + 1;
774 template <
class T,
int N>
776 fftwCorrespondingShapeC2R(TinyVector<T, N> shape,
bool oddDimension0 =
false)
778 shape[0] = oddDimension0
779 ? (shape[0] - 1) * 2 + 1
780 : (shape[0] - 1) * 2;
817 template <
unsigned int N,
class Real =
double>
821 typedef typename FFTWReal2Complex<Real>::plan_type PlanType;
825 Shape shape, instrides, outstrides;
845 template <
class C1,
class C2>
848 int SIGN,
unsigned int planner_flags = FFTW_ESTIMATE)
851 init(in, out, SIGN, planner_flags);
865 template <
class C1,
class C2>
868 unsigned int planner_flags = FFTW_ESTIMATE)
871 init(in, out, planner_flags);
885 template <
class C1,
class C2>
888 unsigned int planner_flags = FFTW_ESTIMATE)
891 init(in, out, planner_flags);
902 instrides.swap(o.instrides);
903 outstrides.swap(o.outstrides);
916 instrides.swap(o.instrides);
917 outstrides.swap(o.outstrides);
928 detail::fftwPlanDestroy(plan);
935 template <
class C1,
class C2>
938 int SIGN,
unsigned int planner_flags = FFTW_ESTIMATE)
940 vigra_precondition(in.strideOrdering() == out.strideOrdering(),
941 "FFTWPlan.init(): input and output must have the same stride ordering.");
943 initImpl(in.permuteStridesDescending(), out.permuteStridesDescending(),
944 SIGN, planner_flags);
951 template <
class C1,
class C2>
954 unsigned int planner_flags = FFTW_ESTIMATE)
957 "FFTWPlan.init(): input and output must have the same stride ordering.");
960 FFTW_FORWARD, planner_flags);
967 template <
class C1,
class C2>
970 unsigned int planner_flags = FFTW_ESTIMATE)
973 "FFTWPlan.init(): input and output must have the same stride ordering.");
976 FFTW_BACKWARD, planner_flags);
986 template <
class C1,
class C2>
990 executeImpl(in.permuteStridesDescending(), out.permuteStridesDescending());
1000 template <
class C1,
class C2>
1014 template <
class C1,
class C2>
1023 template <
class MI,
class MO>
1024 void initImpl(MI ins, MO outs,
int SIGN,
unsigned int planner_flags);
1026 template <
class MI,
class MO>
1027 void executeImpl(MI ins, MO outs)
const;
1032 vigra_precondition(in.shape() == out.shape(),
1033 "FFTWPlan.init(): input and output must have the same shape.");
1036 void checkShapes(MultiArrayView<N, Real, StridedArrayTag> ins,
1037 MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> outs)
const
1039 for(
int k=0; k<(int)N-1; ++k)
1040 vigra_precondition(ins.shape(k) == outs.shape(k),
1041 "FFTWPlan.init(): input and output must have matching shapes.");
1042 vigra_precondition(ins.shape(N-1) / 2 + 1 == outs.shape(N-1),
1043 "FFTWPlan.init(): input and output must have matching shapes.");
1046 void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> ins,
1047 MultiArrayView<N, Real, StridedArrayTag> outs)
const
1049 for(
int k=0; k<(int)N-1; ++k)
1050 vigra_precondition(ins.shape(k) == outs.shape(k),
1051 "FFTWPlan.init(): input and output must have matching shapes.");
1052 vigra_precondition(outs.shape(N-1) / 2 + 1 == ins.shape(N-1),
1053 "FFTWPlan.init(): input and output must have matching shapes.");
1057 template <
unsigned int N,
class Real>
1058 template <
class MI,
class MO>
1060 FFTWPlan<N, Real>::initImpl(MI ins, MO outs,
int SIGN,
unsigned int planner_flags)
1062 checkShapes(ins, outs);
1064 typename MultiArrayShape<N>::type logicalShape(SIGN == FFTW_FORWARD
1068 Shape newShape(logicalShape.begin(), logicalShape.end()),
1069 newIStrides(ins.stride().begin(), ins.stride().end()),
1070 newOStrides(outs.stride().begin(), outs.stride().end()),
1071 itotal(ins.shape().begin(), ins.shape().end()),
1072 ototal(outs.shape().begin(), outs.shape().end());
1074 for(
unsigned int j=1; j<N; ++j)
1076 itotal[j] = ins.stride(j-1) / ins.stride(j);
1077 ototal[j] = outs.stride(j-1) / outs.stride(j);
1080 PlanType newPlan = detail::fftwPlanCreate(N, newShape.begin(),
1081 ins.data(), itotal.begin(), ins.stride(N-1),
1082 outs.data(), ototal.begin(), outs.stride(N-1),
1083 SIGN, planner_flags);
1084 detail::fftwPlanDestroy(plan);
1086 shape.swap(newShape);
1087 instrides.swap(newIStrides);
1088 outstrides.swap(newOStrides);
1092 template <
unsigned int N,
class Real>
1093 template <
class MI,
class MO>
1094 void FFTWPlan<N, Real>::executeImpl(MI ins, MO outs)
const
1096 vigra_precondition(plan != 0,
"FFTWPlan::execute(): plan is NULL.");
1098 typename MultiArrayShape<N>::type lshape(sign == FFTW_FORWARD
1102 vigra_precondition((lshape == TinyVectorView<int, N>(shape.data())),
1103 "FFTWPlan::execute(): shape mismatch between plan and data.");
1104 vigra_precondition((ins.stride() == TinyVectorView<int, N>(instrides.data())),
1105 "FFTWPlan::execute(): strides mismatch between plan and input data.");
1106 vigra_precondition((outs.stride() == TinyVectorView<int, N>(outstrides.data())),
1107 "FFTWPlan::execute(): strides mismatch between plan and output data.");
1109 detail::fftwPlanExecute(plan, ins.data(), outs.data());
1111 typedef typename MO::value_type V;
1112 if(sign == FFTW_BACKWARD)
1113 outs *= V(1.0) / Real(outs.size());
1155 template <
unsigned int N,
class Real>
1163 RArray realArray, realKernel;
1164 CArray fourierArray, fourierKernel;
1165 bool useFourierKernel;
1176 : useFourierKernel(false)
1190 template <
class C1,
class C2,
class C3>
1194 unsigned int planner_flags = FFTW_ESTIMATE)
1195 : useFourierKernel(false)
1197 init(in, kernel, out, planner_flags);
1211 template <
class C1,
class C2,
class C3>
1215 unsigned int planner_flags = FFTW_ESTIMATE)
1216 : useFourierKernel(true)
1218 init(in, kernel, out, planner_flags);
1233 template <
class C1,
class C2,
class C3>
1237 bool fourierDomainKernel,
1238 unsigned int planner_flags = FFTW_ESTIMATE)
1240 init(in, kernel, out, fourierDomainKernel, planner_flags);
1256 template <
class C1,
class C2,
class C3>
1258 bool useFourierKernel =
false,
1259 unsigned int planner_flags = FFTW_ESTIMATE)
1261 if(useFourierKernel)
1262 init(inOut, kernel, planner_flags);
1264 initFourierKernel(inOut, kernel, planner_flags);
1271 template <
class C1,
class C2,
class C3>
1275 unsigned int planner_flags = FFTW_ESTIMATE)
1277 vigra_precondition(in.
shape() == out.
shape(),
1278 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1286 template <
class C1,
class C2,
class C3>
1290 unsigned int planner_flags = FFTW_ESTIMATE)
1292 vigra_precondition(in.
shape() == out.
shape(),
1293 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1294 initFourierKernel(in.
shape(), kernel.shape(), planner_flags);
1301 template <
class C1,
class C2,
class C3>
1305 bool fourierDomainKernel,
1306 unsigned int planner_flags = FFTW_ESTIMATE)
1308 vigra_precondition(in.shape() == out.shape(),
1309 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1310 useFourierKernel = fourierDomainKernel;
1311 initComplex(in.shape(), kernel.shape(), planner_flags);
1320 template <
class C1,
class KernelIterator,
class OutIterator>
1322 KernelIterator kernels, KernelIterator kernelsEnd,
1323 OutIterator outs,
unsigned int planner_flags = FFTW_ESTIMATE)
1325 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1326 typedef typename KernelArray::value_type KernelValue;
1327 typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1328 typedef typename OutArray::value_type OutValue;
1330 bool realKernel = IsSameType<KernelValue, Real>::value;
1331 bool fourierKernel = IsSameType<KernelValue, Complex>::value;
1333 vigra_precondition(realKernel || fourierKernel,
1334 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
1335 vigra_precondition((IsSameType<OutValue, Real>::value),
1336 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
1345 initFourierKernelMany(in.
shape(),
1346 checkShapesFourier(in.
shape(), kernels, kernelsEnd, outs),
1357 template <
class C1,
class KernelIterator,
class OutIterator>
1359 KernelIterator kernels, KernelIterator kernelsEnd,
1361 bool fourierDomainKernels,
1362 unsigned int planner_flags = FFTW_ESTIMATE)
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;
1369 vigra_precondition((IsSameType<KernelValue, Complex>::value),
1370 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
1371 vigra_precondition((IsSameType<OutValue, Complex>::value),
1372 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
1374 useFourierKernel = fourierDomainKernels;
1376 Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs);
1378 CArray newFourierArray(paddedShape), newFourierKernel(paddedShape);
1380 FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags);
1381 FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags);
1383 forward_plan = fplan;
1384 backward_plan = bplan;
1385 fourierArray.
swap(newFourierArray);
1386 fourierKernel.
swap(newFourierKernel);
1389 void init(Shape inOut, Shape kernel,
1390 unsigned int planner_flags = FFTW_ESTIMATE);
1392 void initFourierKernel(Shape inOut, Shape kernel,
1393 unsigned int planner_flags = FFTW_ESTIMATE);
1395 void initComplex(Shape inOut, Shape kernel,
1396 unsigned int planner_flags = FFTW_ESTIMATE);
1398 void initMany(Shape inOut, Shape maxKernel,
1399 unsigned int planner_flags = FFTW_ESTIMATE)
1401 init(inOut, maxKernel, planner_flags);
1404 void initFourierKernelMany(Shape inOut, Shape kernels,
1405 unsigned int planner_flags = FFTW_ESTIMATE)
1407 initFourierKernel(inOut, kernels, planner_flags);
1417 template <
class C1,
class C2,
class C3>
1418 void execute(MultiArrayView<N, Real, C1> in,
1419 MultiArrayView<N, Real, C2> kernel,
1420 MultiArrayView<N, Real, C3> out);
1429 template <
class C1,
class C2,
class C3>
1430 void execute(MultiArrayView<N, Real, C1> in,
1431 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1432 MultiArrayView<N, Real, C3> out);
1441 template <
class C1,
class C2,
class C3>
1442 void execute(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1443 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1444 MultiArrayView<N, FFTWComplex<Real>, C3> out);
1454 template <
class C1,
class KernelIterator,
class OutIterator>
1455 void executeMany(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1456 KernelIterator kernels, KernelIterator kernelsEnd,
1466 template <
class C1,
class KernelIterator,
class OutIterator>
1468 KernelIterator kernels, KernelIterator kernelsEnd,
1471 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1472 typedef typename KernelArray::value_type KernelValue;
1473 typedef typename IsSameType<KernelValue, Complex>::type UseFourierKernel;
1474 typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1475 typedef typename OutArray::value_type OutValue;
1477 bool realKernel = IsSameType<KernelValue, Real>::value;
1478 bool fourierKernel = IsSameType<KernelValue, Complex>::value;
1480 vigra_precondition(realKernel || fourierKernel,
1481 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type.");
1482 vigra_precondition((IsSameType<OutValue, Real>::value),
1483 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type.");
1485 executeManyImpl(in, kernels, kernelsEnd, outs, UseFourierKernel());
1490 template <
class KernelIterator,
class OutIterator>
1491 Shape checkShapes(Shape in,
1492 KernelIterator kernels, KernelIterator kernelsEnd,
1495 template <
class KernelIterator,
class OutIterator>
1496 Shape checkShapesFourier(Shape in,
1497 KernelIterator kernels, KernelIterator kernelsEnd,
1500 template <
class KernelIterator,
class OutIterator>
1501 Shape checkShapesComplex(Shape in,
1502 KernelIterator kernels, KernelIterator kernelsEnd,
1505 template <
class C1,
class KernelIterator,
class OutIterator>
1508 KernelIterator kernels, KernelIterator kernelsEnd,
1509 OutIterator outs, VigraFalseType );
1511 template <
class C1,
class KernelIterator,
class OutIterator>
1514 KernelIterator kernels, KernelIterator kernelsEnd,
1515 OutIterator outs, VigraTrueType );
1519 template <
unsigned int N,
class Real>
1522 unsigned int planner_flags)
1524 Shape paddedShape = fftwBestPaddedShapeR2C(in + kernel - Shape(1)),
1527 CArray newFourierArray(complexShape), newFourierKernel(complexShape);
1529 Shape realStrides = 2*newFourierArray.stride();
1531 RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data());
1532 RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data());
1534 FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags);
1535 FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags);
1537 forward_plan = fplan;
1538 backward_plan = bplan;
1539 realArray = newRealArray;
1540 realKernel = newRealKernel;
1541 fourierArray.swap(newFourierArray);
1542 fourierKernel.swap(newFourierKernel);
1543 useFourierKernel =
false;
1546 template <
unsigned int N,
class Real>
1548 FFTWConvolvePlan<N, Real>::initFourierKernel(Shape in, Shape kernel,
1549 unsigned int planner_flags)
1551 Shape complexShape = kernel,
1552 paddedShape = fftwCorrespondingShapeC2R(complexShape);
1554 for(
unsigned int k=0; k<N; ++k)
1555 vigra_precondition(in[k] <= paddedShape[k],
1556 "FFTWConvolvePlan::init(): kernel too small for given input.");
1558 CArray newFourierArray(complexShape), newFourierKernel(complexShape);
1560 Shape realStrides = 2*newFourierArray.stride();
1562 RArray newRealArray(paddedShape, realStrides, (Real*)newFourierArray.data());
1563 RArray newRealKernel(paddedShape, realStrides, (Real*)newFourierKernel.data());
1565 FFTWPlan<N, Real> fplan(newRealArray, newFourierArray, planner_flags);
1566 FFTWPlan<N, Real> bplan(newFourierArray, newRealArray, planner_flags);
1568 forward_plan = fplan;
1569 backward_plan = bplan;
1570 realArray = newRealArray;
1571 realKernel = newRealKernel;
1572 fourierArray.swap(newFourierArray);
1573 fourierKernel.swap(newFourierKernel);
1574 useFourierKernel =
true;
1577 template <
unsigned int N,
class Real>
1579 FFTWConvolvePlan<N, Real>::initComplex(Shape in, Shape kernel,
1580 unsigned int planner_flags)
1584 if(useFourierKernel)
1586 for(
unsigned int k=0; k<N; ++k)
1587 vigra_precondition(in[k] <= kernel[k],
1588 "FFTWConvolvePlan::init(): kernel too small for given input.");
1590 paddedShape = kernel;
1594 paddedShape = fftwBestPaddedShape(in + kernel - Shape(1));
1597 CArray newFourierArray(paddedShape), newFourierKernel(paddedShape);
1599 FFTWPlan<N, Real> fplan(newFourierArray, newFourierArray, FFTW_FORWARD, planner_flags);
1600 FFTWPlan<N, Real> bplan(newFourierArray, newFourierArray, FFTW_BACKWARD, planner_flags);
1602 forward_plan = fplan;
1603 backward_plan = bplan;
1604 fourierArray.swap(newFourierArray);
1605 fourierKernel.swap(newFourierKernel);
1608 #ifndef DOXYGEN // doxygen documents these functions as free functions
1610 template <
unsigned int N,
class Real>
1611 template <
class C1,
class C2,
class C3>
1614 MultiArrayView<N, Real, C2> kernel,
1615 MultiArrayView<N, Real, C3> out)
1617 vigra_precondition(!useFourierKernel,
1618 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel.");
1620 vigra_precondition(in.shape() == out.shape(),
1621 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1623 Shape paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernel.shape() - Shape(1)),
1624 diff = paddedShape - in.shape(),
1626 right = in.shape() + left;
1628 vigra_precondition(paddedShape == realArray.shape(),
1629 "FFTWConvolvePlan::execute(): shape mismatch between input and plan.");
1631 detail::fftEmbedArray(in, realArray);
1632 forward_plan.execute(realArray, fourierArray);
1634 detail::fftEmbedKernel(kernel, realKernel);
1635 forward_plan.execute(realKernel, fourierKernel);
1637 fourierArray *= fourierKernel;
1639 backward_plan.execute(fourierArray, realArray);
1641 out = realArray.subarray(left, right);
1644 template <
unsigned int N,
class Real>
1645 template <
class C1,
class C2,
class C3>
1648 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1649 MultiArrayView<N, Real, C3> out)
1651 vigra_precondition(useFourierKernel,
1652 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel.");
1654 vigra_precondition(in.shape() == out.shape(),
1655 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1657 vigra_precondition(kernel.shape() == fourierArray.shape(),
1658 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan.");
1660 Shape paddedShape = fftwCorrespondingShapeC2R(kernel.shape(),
odd(in.shape(0))),
1661 diff = paddedShape - in.shape(),
1663 right = in.shape() + left;
1665 vigra_precondition(paddedShape == realArray.shape(),
1666 "FFTWConvolvePlan::execute(): shape mismatch between input and plan.");
1668 detail::fftEmbedArray(in, realArray);
1669 forward_plan.execute(realArray, fourierArray);
1671 fourierKernel = kernel;
1672 moveDCToHalfspaceUpperLeft(fourierKernel);
1674 fourierArray *= fourierKernel;
1676 backward_plan.execute(fourierArray, realArray);
1678 out = realArray.subarray(left, right);
1681 template <
unsigned int N,
class Real>
1682 template <
class C1,
class C2,
class C3>
1685 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1686 MultiArrayView<N, FFTWComplex<Real>, C3> out)
1688 vigra_precondition(in.shape() == out.shape(),
1689 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1691 Shape paddedShape = fourierArray.shape(),
1692 diff = paddedShape - in.shape(),
1694 right = in.shape() + left;
1696 if(useFourierKernel)
1698 vigra_precondition(kernel.shape() == fourierArray.shape(),
1699 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan.");
1701 fourierKernel = kernel;
1706 detail::fftEmbedKernel(kernel, fourierKernel);
1707 forward_plan.execute(fourierKernel, fourierKernel);
1710 detail::fftEmbedArray(in, fourierArray);
1711 forward_plan.execute(fourierArray, fourierArray);
1713 fourierArray *= fourierKernel;
1715 backward_plan.execute(fourierArray, fourierArray);
1717 out = fourierArray.subarray(left, right);
1720 template <
unsigned int N,
class Real>
1721 template <
class C1,
class KernelIterator,
class OutIterator>
1723 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
1724 KernelIterator kernels, KernelIterator kernelsEnd,
1725 OutIterator outs, VigraFalseType )
1727 vigra_precondition(!useFourierKernel,
1728 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel.");
1730 Shape kernelMax = checkShapes(in.shape(), kernels, kernelsEnd, outs),
1731 paddedShape = fftwBestPaddedShapeR2C(in.shape() + kernelMax - Shape(1)),
1732 diff = paddedShape - in.shape(),
1734 right = in.shape() + left;
1736 vigra_precondition(paddedShape == realArray.shape(),
1737 "FFTWConvolvePlan::executeMany(): shape mismatch between input and plan.");
1739 detail::fftEmbedArray(in, realArray);
1740 forward_plan.execute(realArray, fourierArray);
1742 for(; kernels != kernelsEnd; ++kernels, ++outs)
1744 detail::fftEmbedKernel(*kernels, realKernel);
1745 forward_plan.execute(realKernel, fourierKernel);
1747 fourierKernel *= fourierArray;
1749 backward_plan.execute(fourierKernel, realKernel);
1751 *outs = realKernel.subarray(left, right);
1755 template <
unsigned int N,
class Real>
1756 template <
class C1,
class KernelIterator,
class OutIterator>
1758 FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
1759 KernelIterator kernels, KernelIterator kernelsEnd,
1760 OutIterator outs, VigraTrueType )
1762 vigra_precondition(useFourierKernel,
1763 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel.");
1765 Shape complexShape = checkShapesFourier(in.shape(), kernels, kernelsEnd, outs),
1766 paddedShape = fftwCorrespondingShapeC2R(complexShape,
odd(in.shape(0))),
1767 diff = paddedShape - in.shape(),
1769 right = in.shape() + left;
1771 vigra_precondition(complexShape == fourierArray.shape(),
1772 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between kernels and plan.");
1774 vigra_precondition(paddedShape == realArray.shape(),
1775 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between input and plan.");
1777 detail::fftEmbedArray(in, realArray);
1778 forward_plan.execute(realArray, fourierArray);
1780 for(; kernels != kernelsEnd; ++kernels, ++outs)
1782 fourierKernel = *kernels;
1783 moveDCToHalfspaceUpperLeft(fourierKernel);
1784 fourierKernel *= fourierArray;
1786 backward_plan.execute(fourierKernel, realKernel);
1788 *outs = realKernel.subarray(left, right);
1792 template <
unsigned int N,
class Real>
1793 template <
class C1,
class KernelIterator,
class OutIterator>
1796 KernelIterator kernels, KernelIterator kernelsEnd,
1799 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1800 typedef typename KernelArray::value_type KernelValue;
1801 typedef typename std::iterator_traits<OutIterator>::value_type OutArray;
1802 typedef typename OutArray::value_type OutValue;
1804 vigra_precondition((IsSameType<KernelValue, Complex>::value),
1805 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type.");
1806 vigra_precondition((IsSameType<OutValue, Complex>::value),
1807 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type.");
1809 Shape paddedShape = checkShapesComplex(in.shape(), kernels, kernelsEnd, outs),
1810 diff = paddedShape - in.shape(),
1812 right = in.shape() + left;
1814 detail::fftEmbedArray(in, fourierArray);
1815 forward_plan.execute(fourierArray, fourierArray);
1817 for(; kernels != kernelsEnd; ++kernels, ++outs)
1819 if(useFourierKernel)
1821 fourierKernel = *kernels;
1826 detail::fftEmbedKernel(*kernels, fourierKernel);
1827 forward_plan.execute(fourierKernel, fourierKernel);
1830 fourierKernel *= fourierArray;
1832 backward_plan.execute(fourierKernel, fourierKernel);
1834 *outs = fourierKernel.subarray(left, right);
1840 template <
unsigned int N,
class Real>
1841 template <
class KernelIterator,
class OutIterator>
1842 typename FFTWConvolvePlan<N, Real>::Shape
1843 FFTWConvolvePlan<N, Real>::checkShapes(Shape in,
1844 KernelIterator kernels, KernelIterator kernelsEnd,
1847 vigra_precondition(kernels != kernelsEnd,
1848 "FFTWConvolvePlan::checkShapes(): empty kernel sequence.");
1851 for(; kernels != kernelsEnd; ++kernels, ++outs)
1853 vigra_precondition(in == outs->shape(),
1854 "FFTWConvolvePlan::checkShapes(): shape mismatch between input and (one) output.");
1855 kernelMax = max(kernelMax, kernels->shape());
1857 vigra_precondition(
prod(kernelMax) > 0,
1858 "FFTWConvolvePlan::checkShapes(): all kernels have size 0.");
1862 template <
unsigned int N,
class Real>
1863 template <
class KernelIterator,
class OutIterator>
1864 typename FFTWConvolvePlan<N, Real>::Shape
1865 FFTWConvolvePlan<N, Real>::checkShapesFourier(Shape in,
1866 KernelIterator kernels, KernelIterator kernelsEnd,
1869 vigra_precondition(kernels != kernelsEnd,
1870 "FFTWConvolvePlan::checkShapesFourier(): empty kernel sequence.");
1872 Shape complexShape = kernels->shape(),
1873 paddedShape = fftwCorrespondingShapeC2R(complexShape);
1875 for(
unsigned int k=0; k<N; ++k)
1876 vigra_precondition(in[k] <= paddedShape[k],
1877 "FFTWConvolvePlan::checkShapesFourier(): kernels too small for given input.");
1879 for(; kernels != kernelsEnd; ++kernels, ++outs)
1881 vigra_precondition(in == outs->shape(),
1882 "FFTWConvolvePlan::checkShapesFourier(): shape mismatch between input and (one) output.");
1883 vigra_precondition(complexShape == kernels->shape(),
1884 "FFTWConvolvePlan::checkShapesFourier(): all kernels must have the same size.");
1886 return complexShape;
1889 template <
unsigned int N,
class Real>
1890 template <
class KernelIterator,
class OutIterator>
1891 typename FFTWConvolvePlan<N, Real>::Shape
1892 FFTWConvolvePlan<N, Real>::checkShapesComplex(Shape in,
1893 KernelIterator kernels, KernelIterator kernelsEnd,
1896 vigra_precondition(kernels != kernelsEnd,
1897 "FFTWConvolvePlan::checkShapesComplex(): empty kernel sequence.");
1899 Shape kernelShape = kernels->shape();
1900 for(; kernels != kernelsEnd; ++kernels, ++outs)
1902 vigra_precondition(in == outs->shape(),
1903 "FFTWConvolvePlan::checkShapesComplex(): shape mismatch between input and (one) output.");
1904 if(useFourierKernel)
1906 vigra_precondition(kernelShape == kernels->shape(),
1907 "FFTWConvolvePlan::checkShapesComplex(): Fourier domain kernels must have identical size.");
1911 kernelShape = max(kernelShape, kernels->shape());
1914 vigra_precondition(
prod(kernelShape) > 0,
1915 "FFTWConvolvePlan::checkShapesComplex(): all kernels have size 0.");
1917 if(useFourierKernel)
1919 for(
unsigned int k=0; k<N; ++k)
1920 vigra_precondition(in[k] <= kernelShape[k],
1921 "FFTWConvolvePlan::checkShapesComplex(): kernels too small for given input.");
1926 return fftwBestPaddedShape(in + kernelShape - Shape(1));
1937 template <
unsigned int N,
class Real,
class C1,
class C2>
1940 MultiArrayView<N, FFTWComplex<Real>, C2> out)
1942 FFTWPlan<N, Real>(in, out, FFTW_FORWARD).execute(in, out);
1945 template <
unsigned int N,
class Real,
class C1,
class C2>
1948 MultiArrayView<N, FFTWComplex<Real>, C2> out)
1950 FFTWPlan<N, Real>(in, out, FFTW_BACKWARD).execute(in, out);
1953 template <
unsigned int N,
class Real,
class C1,
class C2>
1956 MultiArrayView<N, FFTWComplex<Real>, C2> out)
1958 if(in.shape() == out.shape())
1962 FFTWPlan<N, Real>(out, out, FFTW_FORWARD).execute(out, out);
1966 FFTWPlan<N, Real>(in, out).execute(in, out);
1969 vigra_precondition(
false,
1970 "fourierTransform(): shape mismatch between input and output.");
1973 template <
unsigned int N,
class Real,
class C1,
class C2>
1976 MultiArrayView<N, Real, C2> out)
1979 "fourierTransformInverse(): shape mismatch between input and output.");
1980 FFTWPlan<N, Real>(in, out).execute(in, out);
2172 doxygen_overloaded_function(template <...>
void convolveFFT)
2174 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2177 MultiArrayView<N, Real, C2> kernel,
2178 MultiArrayView<N, Real, C3> out)
2180 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out);
2183 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2186 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2187 MultiArrayView<N, Real, C3> out)
2189 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out);
2198 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2201 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2202 MultiArrayView<N, FFTWComplex<Real>, C3> out,
2203 bool fourierDomainKernel)
2205 FFTWConvolvePlan<N, Real>(in, kernel, out, fourierDomainKernel).execute(in, kernel, out);
2214 template <
unsigned int N,
class Real,
class C1,
2215 class KernelIterator,
class OutIterator>
2218 KernelIterator kernels, KernelIterator kernelsEnd,
2221 FFTWConvolvePlan<N, Real> plan;
2222 plan.initMany(in, kernels, kernelsEnd, outs);
2223 plan.executeMany(in, kernels, kernelsEnd, outs);
2232 template <
unsigned int N,
class Real,
class C1,
2233 class KernelIterator,
class OutIterator>
2236 KernelIterator kernels, KernelIterator kernelsEnd,
2238 bool fourierDomainKernel)
2240 FFTWConvolvePlan<N, Real> plan;
2241 plan.initMany(in, kernels, kernelsEnd, outs, fourierDomainKernel);
2242 plan.executeMany(in, kernels, kernelsEnd, outs);
2249 #endif // VIGRA_MULTI_FFT_HXX