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)
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)
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,
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)
708 Shape diff = out.shape() - in.shape(),
710 rightDiff = diff - leftDiff,
711 right = in.shape() + leftDiff;
713 out.subarray(leftDiff, right) = in;
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.");
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.");
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>
1062 checkShapes(ins, outs);
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>
1096 vigra_precondition(plan != 0,
"FFTWPlan::execute(): plan is NULL.");
1103 "FFTWPlan::execute(): shape mismatch between plan and data.");
1105 "FFTWPlan::execute(): strides mismatch between plan and input 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.");
1279 init(in.
shape(), kernel.
shape(), planner_flags);
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.");
1340 initMany(in.
shape(), checkShapes(in.
shape(), kernels, kernelsEnd, outs),
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>
1429 template <
class C1,
class C2,
class C3>
1441 template <
class C1,
class C2,
class C3>
1454 template <
class C1,
class KernelIterator,
class OutIterator>
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());
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>
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());
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>
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>
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>
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>
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>
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>
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>
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>
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>
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>
1945 template <
unsigned int N,
class Real,
class C1,
class C2>
1953 template <
unsigned int N,
class Real,
class C1,
class C2>
1969 vigra_precondition(
false,
1970 "fourierTransform(): shape mismatch between input and output.");
1973 template <
unsigned int N,
class Real,
class C1,
class C2>
1979 "fourierTransformInverse(): shape mismatch between input and output.");
2172 doxygen_overloaded_function(template <...>
void convolveFFT)
2174 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2183 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2198 template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2203 bool fourierDomainKernel)
2214 template <
unsigned int N,
class Real,
class C1,
2215 class KernelIterator,
class OutIterator>
2218 KernelIterator kernels, KernelIterator kernelsEnd,
2222 plan.
initMany(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)
2241 plan.
initMany(in, kernels, kernelsEnd, outs, fourierDomainKernel);
2249 #endif // VIGRA_MULTI_FFT_HXX void execute(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out) const
Execute a complex-to-complex transform.
Definition: multi_fft.hxx:987
void convolveFFTMany(...)
Convolve a real-valued array with a sequence of kernels by means of the Fourier transform.
void execute(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out)
Execute a plan to convolve a real array with a real kernel.
void convolveFFT(...)
Convolve an array with a kernel by means of the Fourier transform.
~FFTWPlan()
Destructor.
Definition: multi_fft.hxx:926
FFTWReal2Complex< Real >::type complex_type
Definition: fftw3.hxx:136
void init(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, Real, C2 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a complex-to-real transform.
Definition: multi_fft.hxx:968
const difference_type & shape() const
Definition: multi_array.hxx:1551
void init(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a real array with a complex kernel.
Definition: multi_fft.hxx:1287
FFTWConvolvePlan(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, FFTWComplex< Real >, C3 > out, bool fourierDomainKernel, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan to convolve a complex array with a complex kernel.
Definition: multi_fft.hxx:1234
pointer data() const
Definition: multi_array.hxx:1792
void init(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a real array with a real kernel.
Definition: multi_fft.hxx:1272
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2357
void execute(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out) const
Execute a real-to-complex transform.
Definition: multi_fft.hxx:1001
void initMany(MultiArrayView< N, Real, C1 > in, KernelIterator kernels, KernelIterator kernelsEnd, OutIterator outs, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a real array with a sequence of kernels.
Definition: multi_fft.hxx:1321
std::ptrdiff_t MultiArrayIndex
Definition: multi_shape.hxx:55
void init(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out, int SIGN, unsigned int planner_flags=FFTW_ESTIMATE)
Init a complex-to-complex transform.
Definition: multi_fft.hxx:936
void convolveFFTComplexMany(...)
Convolve a complex-valued array with a sequence of kernels by means of the Fourier transform...
Definition: accessor.hxx:43
bool odd(int t)
Check if an integer is odd.
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
vigra::detail::MultiIteratorChooser< StrideTag >::template Traverser< actual_dimension, T, T &, T * >::type traverser
Definition: multi_array.hxx:715
NumericTraits< V >::Promote prod(TinyVectorBase< V, SIZE, D1, D2 > const &l)
product of the vector's elements
Definition: tinyvector.hxx:1895
void executeMany(MultiArrayView< N, FFTWComplex< Real >, C1 > in, KernelIterator kernels, KernelIterator kernelsEnd, OutIterator outs)
Execute a plan to convolve a complex array with a sequence of kernels.
FFTWPlan(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out, int SIGN, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan for a complex-to-complex transform.
Definition: multi_fft.hxx:846
void init(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, FFTWComplex< Real >, C3 > out, bool fourierDomainKernel, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a complex array with a complex kernel.
Definition: multi_fft.hxx:1302
FFTWPlan(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan for a real-to-complex transform.
Definition: multi_fft.hxx:866
Definition: multi_fft.hxx:818
void initMany(MultiArrayView< N, FFTWComplex< Real >, C1 > in, KernelIterator kernels, KernelIterator kernelsEnd, OutIterator outs, bool fourierDomainKernels, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a complex array with a sequence of kernels.
Definition: multi_fft.hxx:1358
MultiArrayView< N, T, StridedArrayTag > permuteStridesDescending() const
Definition: multi_array.hxx:2056
Definition: metaprogramming.hxx:100
TinyVector< MultiArrayIndex, N > type
Definition: multi_shape.hxx:237
bool even(int t)
Check if an integer is even.
Wrapper for fixed size vectors.
Definition: tinyvector.hxx:612
void executeMany(MultiArrayView< N, Real, C1 > in, KernelIterator kernels, KernelIterator kernelsEnd, OutIterator outs)
Execute a plan to convolve a real array with a sequence of kernels.
Definition: multi_fft.hxx:1467
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
FixedPoint16< IntBits3, OverflowHandling > & div(FixedPoint16< IntBits1, OverflowHandling > l, FixedPoint16< IntBits2, OverflowHandling > r, FixedPoint16< IntBits3, OverflowHandling > &result)
division with enforced result type.
Definition: fixedpoint.hxx:1616
FFTWConvolvePlan(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan to convolve a real array with a complex kernel.
Definition: multi_fft.hxx:1212
FFTWPlan(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, Real, C2 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan for a complex-to-real transform.
Definition: multi_fft.hxx:886
Definition: multi_fft.hxx:1156
FFTWPlan(FFTWPlan const &other)
Copy constructor.
Definition: multi_fft.hxx:896
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:655
void convolveFFTComplex(...)
Convolve a complex-valued array by means of the Fourier transform.
FFTWPlan()
Create an empty plan.
Definition: multi_fft.hxx:833
FFTWPlan & operator=(FFTWPlan const &other)
Copy assigment.
Definition: multi_fft.hxx:909
MultiArrayView subarray(difference_type p, difference_type q) const
Definition: multi_array.hxx:1431
FFTWConvolvePlan(Shape inOut, Shape kernel, bool useFourierKernel=false, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan from just the shape information.
Definition: multi_fft.hxx:1257
T sign(T t)
The sign function.
Definition: mathutil.hxx:553
void execute(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, Real, C2 > out) const
Execute a complex-to-real transform.
Definition: multi_fft.hxx:1015
void swap(MultiArray &other)
Definition: multi_array.hxx:2922
void init(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a real-to-complex transform.
Definition: multi_fft.hxx:952
FFTWConvolvePlan()
Create an empty plan.
Definition: multi_fft.hxx:1175
Wrapper class for the FFTW complex types 'fftw_complex'.
Definition: fftw3.hxx:131
difference_type strideOrdering() const
Definition: multi_array.hxx:1520
FFTWConvolvePlan(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan to convolve a real array with a real kernel.
Definition: multi_fft.hxx:1191