36#ifndef VIGRA_MULTI_FFT_HXX
37#define VIGRA_MULTI_FFT_HXX
40#include "metaprogramming.hxx"
41#include "multi_array.hxx"
42#include "multi_math.hxx"
43#include "navigator.hxx"
44#include "copyimage.hxx"
45#include "threading.hxx"
61template <
unsigned int N,
class T,
class C>
62void moveDCToCenterImpl(MultiArrayView<N, T, C> a,
unsigned int startDimension)
65 typedef MultiArrayNavigator<Traverser, N> Navigator;
66 typedef typename Navigator::iterator Iterator;
68 for(
unsigned int d = startDimension; d < N; ++d)
70 Navigator nav(a.traverser_begin(), a.shape(), d);
72 for( ; nav.hasMore(); nav++ )
74 Iterator i = nav.begin();
75 int s = nav.end() - i;
80 for(
int k=0; k<s2; ++k)
82 std::swap(i[k], i[k+s2]);
88 for(
int k=0; k<s2; ++k)
99template <
unsigned int N,
class T,
class C>
100void moveDCToUpperLeftImpl(MultiArrayView<N, T, C> a,
unsigned int startDimension)
103 typedef MultiArrayNavigator<Traverser, N> Navigator;
104 typedef typename Navigator::iterator Iterator;
106 for(
unsigned int d = startDimension; d < N; ++d)
108 Navigator nav(a.traverser_begin(), a.shape(), d);
110 for( ; nav.hasMore(); nav++ )
112 Iterator i = nav.begin();
113 int s = nav.end() - i;
118 for(
int k=0; k<s2; ++k)
120 std::swap(i[k], i[k+s2]);
126 for(
int k=s2; k>0; --k)
145template <
unsigned int N,
class T,
class C>
148 detail::moveDCToCenterImpl(a, 0);
151template <
unsigned int N,
class T,
class C>
154 detail::moveDCToUpperLeftImpl(a, 0);
157template <
unsigned int N,
class T,
class C>
158inline void moveDCToHalfspaceCenter(MultiArrayView<N, T, C> a)
160 detail::moveDCToCenterImpl(a, 1);
163template <
unsigned int N,
class T,
class C>
164inline void moveDCToHalfspaceUpperLeft(MultiArrayView<N, T, C> a)
166 detail::moveDCToUpperLeftImpl(a, 1);
172#ifndef VIGRA_SINGLE_THREADED
174template <
int DUMMY=0>
178 threading::lock_guard<threading::mutex> guard_;
181 : guard_(plan_mutex_)
184 static threading::mutex plan_mutex_;
188threading::mutex FFTWLock<DUMMY>::plan_mutex_;
192template <
int DUMMY=0>
204fftwPlanCreate(
unsigned int N,
int* shape,
205 FFTWComplex<double> * in,
int* instrides,
int instep,
206 FFTWComplex<double> * out,
int* outstrides,
int outstep,
207 int sign,
unsigned int planner_flags)
209 return fftw_plan_many_dft(N, shape, 1,
210 (fftw_complex *)in, instrides, instep, 0,
211 (fftw_complex *)out, outstrides, outstep, 0,
212 sign, planner_flags);
216fftwPlanCreate(
unsigned int N,
int* shape,
217 double * in,
int* instrides,
int instep,
218 FFTWComplex<double> * out,
int* outstrides,
int outstep,
219 int ,
unsigned int planner_flags)
221 return fftw_plan_many_dft_r2c(N, shape, 1,
222 in, instrides, instep, 0,
223 (fftw_complex *)out, outstrides, outstep, 0,
228fftwPlanCreate(
unsigned int N,
int* shape,
229 FFTWComplex<double> * in,
int* instrides,
int instep,
230 double * out,
int* outstrides,
int outstep,
231 int ,
unsigned int planner_flags)
233 return fftw_plan_many_dft_c2r(N, shape, 1,
234 (fftw_complex *)in, instrides, instep, 0,
235 out, outstrides, outstep, 0,
240fftwPlanCreate(
unsigned int N,
int* shape,
241 FFTWComplex<float> * in,
int* instrides,
int instep,
242 FFTWComplex<float> * out,
int* outstrides,
int outstep,
243 int sign,
unsigned int planner_flags)
245 return fftwf_plan_many_dft(N, shape, 1,
246 (fftwf_complex *)in, instrides, instep, 0,
247 (fftwf_complex *)out, outstrides, outstep, 0,
248 sign, planner_flags);
252fftwPlanCreate(
unsigned int N,
int* shape,
253 float * in,
int* instrides,
int instep,
254 FFTWComplex<float> * out,
int* outstrides,
int outstep,
255 int ,
unsigned int planner_flags)
257 return fftwf_plan_many_dft_r2c(N, shape, 1,
258 in, instrides, instep, 0,
259 (fftwf_complex *)out, outstrides, outstep, 0,
264fftwPlanCreate(
unsigned int N,
int* shape,
265 FFTWComplex<float> * in,
int* instrides,
int instep,
266 float * out,
int* outstrides,
int outstep,
267 int ,
unsigned int planner_flags)
269 return fftwf_plan_many_dft_c2r(N, shape, 1,
270 (fftwf_complex *)in, instrides, instep, 0,
271 out, outstrides, outstep, 0,
276fftwPlanCreate(
unsigned int N,
int* shape,
277 FFTWComplex<long double> * in,
int* instrides,
int instep,
278 FFTWComplex<long double> * out,
int* outstrides,
int outstep,
279 int sign,
unsigned int planner_flags)
281 return fftwl_plan_many_dft(N, shape, 1,
282 (fftwl_complex *)in, instrides, instep, 0,
283 (fftwl_complex *)out, outstrides, outstep, 0,
284 sign, planner_flags);
288fftwPlanCreate(
unsigned int N,
int* shape,
289 long double * in,
int* instrides,
int instep,
290 FFTWComplex<long double> * out,
int* outstrides,
int outstep,
291 int ,
unsigned int planner_flags)
293 return fftwl_plan_many_dft_r2c(N, shape, 1,
294 in, instrides, instep, 0,
295 (fftwl_complex *)out, outstrides, outstep, 0,
300fftwPlanCreate(
unsigned int N,
int* shape,
301 FFTWComplex<long double> * in,
int* instrides,
int instep,
302 long double * out,
int* outstrides,
int outstep,
303 int ,
unsigned int planner_flags)
305 return fftwl_plan_many_dft_c2r(N, shape, 1,
306 (fftwl_complex *)in, instrides, instep, 0,
307 out, outstrides, outstep, 0,
311inline void fftwPlanDestroy(fftw_plan plan)
314 fftw_destroy_plan(plan);
317inline void fftwPlanDestroy(fftwf_plan plan)
320 fftwf_destroy_plan(plan);
323inline void fftwPlanDestroy(fftwl_plan plan)
326 fftwl_destroy_plan(plan);
330fftwPlanExecute(fftw_plan plan)
336fftwPlanExecute(fftwf_plan plan)
342fftwPlanExecute(fftwl_plan plan)
348fftwPlanExecute(fftw_plan plan, FFTWComplex<double> * in, FFTWComplex<double> * out)
350 fftw_execute_dft(plan, (fftw_complex *)in, (fftw_complex *)out);
354fftwPlanExecute(fftw_plan plan,
double * in, FFTWComplex<double> * out)
356 fftw_execute_dft_r2c(plan, in, (fftw_complex *)out);
360fftwPlanExecute(fftw_plan plan, FFTWComplex<double> * in,
double * out)
362 fftw_execute_dft_c2r(plan, (fftw_complex *)in, out);
366fftwPlanExecute(fftwf_plan plan, FFTWComplex<float> * in, FFTWComplex<float> * out)
368 fftwf_execute_dft(plan, (fftwf_complex *)in, (fftwf_complex *)out);
372fftwPlanExecute(fftwf_plan plan,
float * in, FFTWComplex<float> * out)
374 fftwf_execute_dft_r2c(plan, in, (fftwf_complex *)out);
378fftwPlanExecute(fftwf_plan plan, FFTWComplex<float> * in,
float * out)
380 fftwf_execute_dft_c2r(plan, (fftwf_complex *)in, out);
384fftwPlanExecute(fftwl_plan plan, FFTWComplex<long double> * in, FFTWComplex<long double> * out)
386 fftwl_execute_dft(plan, (fftwl_complex *)in, (fftwl_complex *)out);
390fftwPlanExecute(fftwl_plan plan,
long double * in, FFTWComplex<long double> * out)
392 fftwl_execute_dft_r2c(plan, in, (fftwl_complex *)out);
396fftwPlanExecute(fftwl_plan plan, FFTWComplex<long double> * in,
long double * out)
398 fftwl_execute_dft_c2r(plan, (fftwl_complex *)in, out);
402struct FFTWPaddingSize
404 static const int size = 1330;
405 static const int evenSize = 1063;
406 static int goodSizes[size];
407 static int goodEvenSizes[evenSize];
409 static inline int find(
int s)
411 if(s <= 0 || s >= goodSizes[size-1])
414 int * upperBound = std::upper_bound(goodSizes, goodSizes+size, s);
415 if(upperBound > goodSizes && upperBound[-1] == s)
421 static inline int findEven(
int s)
423 if(s <= 0 || s >= goodEvenSizes[evenSize-1])
426 int * upperBound = std::upper_bound(goodEvenSizes, goodEvenSizes+evenSize, s);
427 if(upperBound > goodEvenSizes && upperBound[-1] == s)
439int FFTWPaddingSize<DUMMY>::goodSizes[size] = {
440 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
441 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32, 33, 35, 36, 39, 40, 42, 44, 45, 48,
442 49, 50, 52, 54, 55, 56, 60, 63, 64, 65, 66, 70, 72, 75, 77, 78, 80, 81,
443 84, 88, 90, 91, 96, 98, 99, 100, 104, 105, 108, 110, 112, 117, 120, 125,
444 126, 128, 130, 132, 135, 140, 144, 147, 150, 154, 156, 160, 162, 165,
445 168, 175, 176, 180, 182, 189, 192, 195, 196, 198, 200, 208, 210, 216,
446 220, 224, 225, 231, 234, 240, 243, 245, 250, 252, 256, 260, 264, 270,
447 273, 275, 280, 288, 294, 297, 300, 308, 312, 315, 320, 324, 325, 330,
448 336, 343, 350, 351, 352, 360, 364, 375, 378, 384, 385, 390, 392, 396,
449 400, 405, 416, 420, 432, 440, 441, 448, 450, 455, 462, 468, 480, 486,
450 490, 495, 500, 504, 512, 520, 525, 528, 539, 540, 546, 550, 560, 567,
451 576, 585, 588, 594, 600, 616, 624, 625, 630, 637, 640, 648, 650, 660,
452 672, 675, 686, 693, 700, 702, 704, 720, 728, 729, 735, 750, 756, 768,
453 770, 780, 784, 792, 800, 810, 819, 825, 832, 840, 864, 875, 880, 882,
454 891, 896, 900, 910, 924, 936, 945, 960, 972, 975, 980, 990, 1000, 1008,
455 1024, 1029, 1040, 1050, 1053, 1056, 1078, 1080, 1092, 1100, 1120, 1125,
456 1134, 1152, 1155, 1170, 1176, 1188, 1200, 1215, 1225, 1232, 1248, 1250,
457 1260, 1274, 1280, 1296, 1300, 1320, 1323, 1344, 1350, 1365, 1372, 1375,
458 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470, 1485, 1500, 1512, 1536,
459 1540, 1560, 1568, 1575, 1584, 1600, 1617, 1620, 1625, 1638, 1650, 1664,
460 1680, 1701, 1715, 1728, 1750, 1755, 1760, 1764, 1782, 1792, 1800, 1820,
461 1848, 1872, 1875, 1890, 1911, 1920, 1925, 1944, 1950, 1960, 1980, 2000,
462 2016, 2025, 2048, 2058, 2079, 2080, 2100, 2106, 2112, 2156, 2160, 2184,
463 2187, 2200, 2205, 2240, 2250, 2268, 2275, 2304, 2310, 2340, 2352, 2376,
464 2400, 2401, 2430, 2450, 2457, 2464, 2475, 2496, 2500, 2520, 2548, 2560,
465 2592, 2600, 2625, 2640, 2646, 2673, 2688, 2695, 2700, 2730, 2744, 2750,
466 2772, 2800, 2808, 2816, 2835, 2880, 2912, 2916, 2925, 2940, 2970, 3000,
467 3024, 3072, 3080, 3087, 3120, 3125, 3136, 3150, 3159, 3168, 3185, 3200,
468 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3375, 3402, 3430, 3456, 3465,
469 3500, 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3645, 3675, 3696, 3744,
470 3750, 3773, 3780, 3822, 3840, 3850, 3888, 3900, 3920, 3960, 3969, 4000,
471 4032, 4050, 4095, 4096, 4116, 4125, 4158, 4160, 4200, 4212, 4224, 4312,
472 4320, 4368, 4374, 4375, 4400, 4410, 4455, 4459, 4480, 4500, 4536, 4550,
473 4608, 4620, 4680, 4704, 4725, 4752, 4800, 4802, 4851, 4860, 4875, 4900,
474 4914, 4928, 4950, 4992, 5000, 5040, 5096, 5103, 5120, 5145, 5184, 5200,
475 5250, 5265, 5280, 5292, 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544,
476 5600, 5616, 5625, 5632, 5670, 5733, 5760, 5775, 5824, 5832, 5850, 5880,
477 5940, 6000, 6048, 6075, 6125, 6144, 6160, 6174, 6237, 6240, 6250, 6272,
478 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6561, 6600, 6615,
479 6656, 6720, 6750, 6804, 6825, 6860, 6875, 6912, 6930, 7000, 7020, 7040,
480 7056, 7128, 7168, 7200, 7203, 7280, 7290, 7350, 7371, 7392, 7425, 7488,
481 7500, 7546, 7560, 7644, 7680, 7700, 7776, 7800, 7840, 7875, 7920, 7938,
482 8000, 8019, 8064, 8085, 8100, 8125, 8190, 8192, 8232, 8250, 8316, 8320,
483 8400, 8424, 8448, 8505, 8575, 8624, 8640, 8736, 8748, 8750, 8775, 8800,
484 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9261, 9360, 9375,
485 9408, 9450, 9477, 9504, 9555, 9600, 9604, 9625, 9702, 9720, 9750, 9800,
486 9828, 9856, 9900, 9984, 10000, 10080, 10125, 10192, 10206, 10240, 10290,
487 10368, 10395, 10400, 10500, 10530, 10560, 10584, 10692, 10752, 10780,
488 10800, 10920, 10935, 10976, 11000, 11025, 11088, 11200, 11232, 11250,
489 11264, 11319, 11340, 11375, 11466, 11520, 11550, 11648, 11664, 11700,
490 11760, 11880, 11907, 12000, 12005, 12096, 12150, 12250, 12285, 12288,
491 12320, 12348, 12375, 12474, 12480, 12500, 12544, 12600, 12636, 12672,
492 12740, 12800, 12936, 12960, 13000, 13104, 13122, 13125, 13200, 13230,
493 13312, 13365, 13377, 13440, 13475, 13500, 13608, 13650, 13720, 13750,
494 13824, 13860, 14000, 14040, 14080, 14112, 14175, 14256, 14336, 14400,
495 14406, 14553, 14560, 14580, 14625, 14700, 14742, 14784, 14850, 14976,
496 15000, 15092, 15120, 15288, 15309, 15360, 15400, 15435, 15552, 15600,
497 15625, 15680, 15750, 15795, 15840, 15876, 15925, 16000, 16038, 16128,
498 16170, 16200, 16250, 16380, 16384, 16464, 16500, 16632, 16640, 16800,
499 16807, 16848, 16875, 16896, 17010, 17150, 17199, 17248, 17280, 17325,
500 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000,
501 18144, 18200, 18225, 18375, 18432, 18480, 18522, 18711, 18720, 18750,
502 18816, 18865, 18900, 18954, 19008, 19110, 19200, 19208, 19250, 19404,
503 19440, 19500, 19600, 19656, 19683, 19712, 19800, 19845, 19968, 20000,
504 20160, 20250, 20384, 20412, 20475, 20480, 20580, 20625, 20736, 20790,
505 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560, 21600, 21609,
506 21840, 21870, 21875, 21952, 22000, 22050, 22113, 22176, 22275, 22295,
507 22400, 22464, 22500, 22528, 22638, 22680, 22750, 22932, 23040, 23100,
508 23296, 23328, 23400, 23520, 23625, 23760, 23814, 24000, 24010, 24057,
509 24192, 24255, 24300, 24375, 24500, 24570, 24576, 24640, 24696, 24750,
510 24948, 24960, 25000, 25088, 25200, 25272, 25344, 25480, 25515, 25600,
511 25725, 25872, 25920, 26000, 26208, 26244, 26250, 26325, 26400, 26411,
512 26460, 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440,
513 27500, 27648, 27720, 27783, 28000, 28080, 28125, 28160, 28224, 28350,
514 28431, 28512, 28665, 28672, 28800, 28812, 28875, 29106, 29120, 29160,
515 29250, 29400, 29484, 29568, 29700, 29952, 30000, 30184, 30240, 30375,
516 30576, 30618, 30625, 30720, 30800, 30870, 31104, 31185, 31200, 31213,
517 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256,
518 32340, 32400, 32500, 32760, 32768, 32805, 32928, 33000, 33075, 33264,
519 33280, 33600, 33614, 33696, 33750, 33792, 33957, 34020, 34125, 34300,
520 34375, 34398, 34496, 34560, 34650, 34944, 34992, 35000, 35100, 35200,
521 35280, 35640, 35672, 35721, 35840, 36000, 36015, 36288, 36400, 36450,
522 36750, 36855, 36864, 36960, 37044, 37125, 37422, 37440, 37500, 37632,
523 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808, 38880,
524 39000, 39200, 39312, 39366, 39375, 39424, 39600, 39690, 39936, 40000,
525 40095, 40131, 40320, 40425, 40500, 40625, 40768, 40824, 40950, 40960,
526 41160, 41250, 41472, 41580, 41600, 42000, 42120, 42240, 42336, 42525,
527 42768, 42875, 43008, 43120, 43200, 43218, 43659, 43680, 43740, 43750,
528 43875, 43904, 44000, 44100, 44226, 44352, 44550, 44590, 44800, 44928,
529 45000, 45056, 45276, 45360, 45500, 45864, 45927, 46080, 46200, 46305,
530 46592, 46656, 46800, 46875, 47040, 47250, 47385, 47520, 47628, 47775,
531 48000, 48020, 48114, 48125, 48384, 48510, 48600, 48750, 49000, 49140,
532 49152, 49280, 49392, 49500, 49896, 49920, 50000, 50176, 50400, 50421,
533 50544, 50625, 50688, 50960, 51030, 51200, 51450, 51597, 51744, 51840,
534 51975, 52000, 52416, 52488, 52500, 52650, 52800, 52822, 52920, 53248,
535 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54675, 54880, 55000,
536 55125, 55296, 55440, 55566, 56000, 56133, 56160, 56250, 56320, 56448,
537 56595, 56700, 56862, 56875, 57024, 57330, 57344, 57600, 57624, 57750,
538 58212, 58240, 58320, 58500, 58800, 58968, 59049, 59136, 59400, 59535,
539 59904, 60000, 60025, 60368, 60480, 60750, 61152, 61236, 61250, 61425,
540 61440, 61600, 61740, 61875, 62208, 62370, 62400, 62426, 62500, 62720,
541 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512, 64680, 64800,
542 64827, 65000, 65520, 65536, 65610, 65625, 65856, 66000, 66150, 66339,
543 66528, 66560, 66825, 66885, 67200, 67228, 67375, 67392, 67500, 67584,
544 67914, 68040, 68250, 68600, 68750, 68796, 68992, 69120, 69300, 69888,
545 69984, 70000, 70200, 70400, 70560, 70875, 71280, 71344, 71442, 71680,
546 72000, 72030, 72171, 72576, 72765, 72800, 72900, 73125, 73500, 73710,
547 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264, 75460, 75600,
548 75816, 76032, 76440, 76545, 76800, 76832, 77000, 77175, 77616, 77760,
549 78000, 78125, 78400, 78624, 78732, 78750, 78848, 78975, 79200, 79233,
550 79380, 79625, 79872, 80000, 80190, 80262, 80640, 80850, 81000, 81250,
551 81536, 81648, 81900, 81920, 82320, 82500, 82944, 83160, 83200, 83349,
552 84000, 84035, 84240, 84375, 84480, 84672, 85050, 85293, 85536, 85750,
553 85995, 86016, 86240, 86400, 86436, 86625, 87318, 87360, 87480, 87500,
554 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856,
555 90000, 90112, 90552, 90720, 91000, 91125, 91728, 91854, 91875, 92160,
556 92400, 92610, 93184, 93312, 93555, 93600, 93639, 93750, 94080, 94325,
557 94500, 94770, 95040, 95256, 95550, 96000, 96040, 96228, 96250, 96768,
558 97020, 97200, 97500, 98000, 98280, 98304, 98415, 98560, 98784, 99000,
563int FFTWPaddingSize<DUMMY>::goodEvenSizes[evenSize] = {
564 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22,
565 24, 26, 28, 30, 32, 36, 40, 42, 44, 48, 50, 52, 54, 56, 60, 64, 66, 70,
566 72, 78, 80, 84, 88, 90, 96, 98, 100, 104, 108, 110, 112, 120, 126, 128,
567 130, 132, 140, 144, 150, 154, 156, 160, 162, 168, 176, 180, 182, 192,
568 196, 198, 200, 208, 210, 216, 220, 224, 234, 240, 250, 252, 256, 260,
569 264, 270, 280, 288, 294, 300, 308, 312, 320, 324, 330, 336, 350, 352,
570 360, 364, 378, 384, 390, 392, 396, 400, 416, 420, 432, 440, 448, 450,
571 462, 468, 480, 486, 490, 500, 504, 512, 520, 528, 540, 546, 550, 560,
572 576, 588, 594, 600, 616, 624, 630, 640, 648, 650, 660, 672, 686, 700,
573 702, 704, 720, 728, 750, 756, 768, 770, 780, 784, 792, 800, 810, 832,
574 840, 864, 880, 882, 896, 900, 910, 924, 936, 960, 972, 980, 990, 1000,
575 1008, 1024, 1040, 1050, 1056, 1078, 1080, 1092, 1100, 1120, 1134, 1152,
576 1170, 1176, 1188, 1200, 1232, 1248, 1250, 1260, 1274, 1280, 1296, 1300,
577 1320, 1344, 1350, 1372, 1386, 1400, 1404, 1408, 1440, 1456, 1458, 1470,
578 1500, 1512, 1536, 1540, 1560, 1568, 1584, 1600, 1620, 1638, 1650, 1664,
579 1680, 1728, 1750, 1760, 1764, 1782, 1792, 1800, 1820, 1848, 1872, 1890,
580 1920, 1944, 1950, 1960, 1980, 2000, 2016, 2048, 2058, 2080, 2100, 2106,
581 2112, 2156, 2160, 2184, 2200, 2240, 2250, 2268, 2304, 2310, 2340, 2352,
582 2376, 2400, 2430, 2450, 2464, 2496, 2500, 2520, 2548, 2560, 2592, 2600,
583 2640, 2646, 2688, 2700, 2730, 2744, 2750, 2772, 2800, 2808, 2816, 2880,
584 2912, 2916, 2940, 2970, 3000, 3024, 3072, 3080, 3120, 3136, 3150, 3168,
585 3200, 3234, 3240, 3250, 3276, 3300, 3328, 3360, 3402, 3430, 3456, 3500,
586 3510, 3520, 3528, 3564, 3584, 3600, 3640, 3696, 3744, 3750, 3780, 3822,
587 3840, 3850, 3888, 3900, 3920, 3960, 4000, 4032, 4050, 4096, 4116, 4158,
588 4160, 4200, 4212, 4224, 4312, 4320, 4368, 4374, 4400, 4410, 4480, 4500,
589 4536, 4550, 4608, 4620, 4680, 4704, 4752, 4800, 4802, 4860, 4900, 4914,
590 4928, 4950, 4992, 5000, 5040, 5096, 5120, 5184, 5200, 5250, 5280, 5292,
591 5346, 5376, 5390, 5400, 5460, 5488, 5500, 5544, 5600, 5616, 5632, 5670,
592 5760, 5824, 5832, 5850, 5880, 5940, 6000, 6048, 6144, 6160, 6174, 6240,
593 6250, 6272, 6300, 6318, 6336, 6370, 6400, 6468, 6480, 6500, 6552, 6600,
594 6656, 6720, 6750, 6804, 6860, 6912, 6930, 7000, 7020, 7040, 7056, 7128,
595 7168, 7200, 7280, 7290, 7350, 7392, 7488, 7500, 7546, 7560, 7644, 7680,
596 7700, 7776, 7800, 7840, 7920, 7938, 8000, 8064, 8100, 8190, 8192, 8232,
597 8250, 8316, 8320, 8400, 8424, 8448, 8624, 8640, 8736, 8748, 8750, 8800,
598 8820, 8910, 8918, 8960, 9000, 9072, 9100, 9216, 9240, 9360, 9408, 9450,
599 9504, 9600, 9604, 9702, 9720, 9750, 9800, 9828, 9856, 9900, 9984, 10000,
600 10080, 10192, 10206, 10240, 10290, 10368, 10400, 10500, 10530, 10560,
601 10584, 10692, 10752, 10780, 10800, 10920, 10976, 11000, 11088, 11200,
602 11232, 11250, 11264, 11340, 11466, 11520, 11550, 11648, 11664, 11700,
603 11760, 11880, 12000, 12096, 12150, 12250, 12288, 12320, 12348, 12474,
604 12480, 12500, 12544, 12600, 12636, 12672, 12740, 12800, 12936, 12960,
605 13000, 13104, 13122, 13200, 13230, 13312, 13440, 13500, 13608, 13650,
606 13720, 13750, 13824, 13860, 14000, 14040, 14080, 14112, 14256, 14336,
607 14400, 14406, 14560, 14580, 14700, 14742, 14784, 14850, 14976, 15000,
608 15092, 15120, 15288, 15360, 15400, 15552, 15600, 15680, 15750, 15840,
609 15876, 16000, 16038, 16128, 16170, 16200, 16250, 16380, 16384, 16464,
610 16500, 16632, 16640, 16800, 16848, 16896, 17010, 17150, 17248, 17280,
611 17472, 17496, 17500, 17550, 17600, 17640, 17820, 17836, 17920, 18000,
612 18144, 18200, 18432, 18480, 18522, 18720, 18750, 18816, 18900, 18954,
613 19008, 19110, 19200, 19208, 19250, 19404, 19440, 19500, 19600, 19656,
614 19712, 19800, 19968, 20000, 20160, 20250, 20384, 20412, 20480, 20580,
615 20736, 20790, 20800, 21000, 21060, 21120, 21168, 21384, 21504, 21560,
616 21600, 21840, 21870, 21952, 22000, 22050, 22176, 22400, 22464, 22500,
617 22528, 22638, 22680, 22750, 22932, 23040, 23100, 23296, 23328, 23400,
618 23520, 23760, 23814, 24000, 24010, 24192, 24300, 24500, 24570, 24576,
619 24640, 24696, 24750, 24948, 24960, 25000, 25088, 25200, 25272, 25344,
620 25480, 25600, 25872, 25920, 26000, 26208, 26244, 26250, 26400, 26460,
621 26624, 26730, 26754, 26880, 26950, 27000, 27216, 27300, 27440, 27500,
622 27648, 27720, 28000, 28080, 28160, 28224, 28350, 28512, 28672, 28800,
623 28812, 29106, 29120, 29160, 29250, 29400, 29484, 29568, 29700, 29952,
624 30000, 30184, 30240, 30576, 30618, 30720, 30800, 30870, 31104, 31200,
625 31250, 31360, 31500, 31590, 31680, 31752, 31850, 32000, 32076, 32256,
626 32340, 32400, 32500, 32760, 32768, 32928, 33000, 33264, 33280, 33600,
627 33614, 33696, 33750, 33792, 34020, 34300, 34398, 34496, 34560, 34650,
628 34944, 34992, 35000, 35100, 35200, 35280, 35640, 35672, 35840, 36000,
629 36288, 36400, 36450, 36750, 36864, 36960, 37044, 37422, 37440, 37500,
630 37632, 37730, 37800, 37908, 38016, 38220, 38400, 38416, 38500, 38808,
631 38880, 39000, 39200, 39312, 39366, 39424, 39600, 39690, 39936, 40000,
632 40320, 40500, 40768, 40824, 40950, 40960, 41160, 41250, 41472, 41580,
633 41600, 42000, 42120, 42240, 42336, 42768, 43008, 43120, 43200, 43218,
634 43680, 43740, 43750, 43904, 44000, 44100, 44226, 44352, 44550, 44590,
635 44800, 44928, 45000, 45056, 45276, 45360, 45500, 45864, 46080, 46200,
636 46592, 46656, 46800, 47040, 47250, 47520, 47628, 48000, 48020, 48114,
637 48384, 48510, 48600, 48750, 49000, 49140, 49152, 49280, 49392, 49500,
638 49896, 49920, 50000, 50176, 50400, 50544, 50688, 50960, 51030, 51200,
639 51450, 51744, 51840, 52000, 52416, 52488, 52500, 52650, 52800, 52822,
640 52920, 53248, 53460, 53508, 53760, 53900, 54000, 54432, 54600, 54880,
641 55000, 55296, 55440, 55566, 56000, 56160, 56250, 56320, 56448, 56700,
642 56862, 57024, 57330, 57344, 57600, 57624, 57750, 58212, 58240, 58320,
643 58500, 58800, 58968, 59136, 59400, 59904, 60000, 60368, 60480, 60750,
644 61152, 61236, 61250, 61440, 61600, 61740, 62208, 62370, 62400, 62426,
645 62500, 62720, 63000, 63180, 63360, 63504, 63700, 64000, 64152, 64512,
646 64680, 64800, 65000, 65520, 65536, 65610, 65856, 66000, 66150, 66528,
647 66560, 67200, 67228, 67392, 67500, 67584, 67914, 68040, 68250, 68600,
648 68750, 68796, 68992, 69120, 69300, 69888, 69984, 70000, 70200, 70400,
649 70560, 71280, 71344, 71442, 71680, 72000, 72030, 72576, 72800, 72900,
650 73500, 73710, 73728, 73920, 74088, 74250, 74844, 74880, 75000, 75264,
651 75460, 75600, 75816, 76032, 76440, 76800, 76832, 77000, 77616, 77760,
652 78000, 78400, 78624, 78732, 78750, 78848, 79200, 79380, 79872, 80000,
653 80190, 80262, 80640, 80850, 81000, 81250, 81536, 81648, 81900, 81920,
654 82320, 82500, 82944, 83160, 83200, 84000, 84240, 84480, 84672, 85050,
655 85536, 85750, 86016, 86240, 86400, 86436, 87318, 87360, 87480, 87500,
656 87750, 87808, 88000, 88200, 88452, 88704, 89100, 89180, 89600, 89856,
657 90000, 90112, 90552, 90720, 91000, 91728, 91854, 92160, 92400, 92610,
658 93184, 93312, 93600, 93750, 94080, 94500, 94770, 95040, 95256, 95550,
659 96000, 96040, 96228, 96250, 96768, 97020, 97200, 97500, 98000, 98280,
660 98304, 98560, 98784, 99000, 99792, 99840
666 template <
unsigned int N,
class Real,
class C,
class Shape>
668 exec(MultiArrayView<N, Real, C> & out, Shape
const & kernelShape,
669 Shape & srcPoint, Shape & destPoint,
bool copyIt)
671 for(srcPoint[M]=0; srcPoint[M]<kernelShape[M]; ++srcPoint[M])
673 if(srcPoint[M] < (kernelShape[M] + 1) / 2)
675 destPoint[M] = srcPoint[M];
679 destPoint[M] = srcPoint[M] + out.shape(M) - kernelShape[M];
682 FFTEmbedKernel<M-1>::exec(out, kernelShape, srcPoint, destPoint, copyIt);
688struct FFTEmbedKernel<0>
690 template <
unsigned int N,
class Real,
class C,
class Shape>
692 exec(MultiArrayView<N, Real, C> & out, Shape
const & kernelShape,
693 Shape & srcPoint, Shape & destPoint,
bool copyIt)
695 for(srcPoint[0]=0; srcPoint[0]<kernelShape[0]; ++srcPoint[0])
697 if(srcPoint[0] < (kernelShape[0] + 1) / 2)
699 destPoint[0] = srcPoint[0];
703 destPoint[0] = srcPoint[0] + out.shape(0) - kernelShape[0];
708 out[destPoint] = out[srcPoint];
715template <
unsigned int N,
class Real,
class C1,
class C2>
717fftEmbedKernel(MultiArrayView<N, Real, C1> kernel,
718 MultiArrayView<N, Real, C2> out,
723 MultiArrayView<N, Real, C2> kout = out.
subarray(Shape(), kernel.shape());
731 Shape srcPoint, destPoint;
732 FFTEmbedKernel<(int)N-1>::exec(out, kernel.shape(), srcPoint, destPoint,
false);
735template <
unsigned int N,
class Real,
class C1,
class C2>
737fftEmbedArray(MultiArrayView<N, Real, C1> in,
738 MultiArrayView<N, Real, C2> out)
742 Shape diff = out.shape() - in.shape(),
744 rightDiff = diff - leftDiff,
745 right = in.shape() + leftDiff;
750 typedef MultiArrayNavigator<Traverser, N> Navigator;
751 typedef typename Navigator::iterator Iterator;
753 for(
unsigned int d = 0; d < N; ++d)
755 Navigator nav(out.traverser_begin(), out.shape(), d);
757 for( ; nav.hasMore(); nav++ )
759 Iterator i = nav.begin();
760 for(
int k=1; k<=leftDiff[d]; ++k)
761 i[leftDiff[d] - k] = i[leftDiff[d] + k];
762 for(
int k=0; k<rightDiff[d]; ++k)
763 i[right[d] + k] = i[right[d] - k - 2];
770template <
class T,
int N>
772fftwBestPaddedShape(TinyVector<T, N> shape)
774 for(
unsigned int k=0; k<N; ++k)
775 shape[k] = detail::FFTWPaddingSize<0>::find(shape[k]);
779template <
class T,
int N>
781fftwBestPaddedShapeR2C(TinyVector<T, N> shape)
783 shape[0] = detail::FFTWPaddingSize<0>::findEven(shape[0]);
784 for(
unsigned int k=1; k<N; ++k)
785 shape[k] = detail::FFTWPaddingSize<0>::find(shape[k]);
800template <
class T,
int N>
804 shape[0] = shape[0] / 2 + 1;
808template <
class T,
int N>
810fftwCorrespondingShapeC2R(TinyVector<T, N> shape,
bool oddDimension0 =
false)
812 shape[0] = oddDimension0
813 ? (shape[0] - 1) * 2 + 1
814 : (shape[0] - 1) * 2;
851template <
unsigned int N,
class Real =
double>
855 typedef typename FFTWReal2Complex<Real>::plan_type
PlanType;
859 Shape shape, instrides, outstrides;
879 template <
class C1,
class C2>
899 template <
class C1,
class C2>
919 template <
class C1,
class C2>
936 instrides.swap(o.instrides);
937 outstrides.swap(o.outstrides);
950 instrides.swap(o.instrides);
951 outstrides.swap(o.outstrides);
962 detail::FFTWLock<> lock;
963 detail::fftwPlanDestroy(plan);
970 template <
class C1,
class C2>
975 vigra_precondition(
in.strideOrdering() ==
out.strideOrdering(),
976 "FFTWPlan.init(): input and output must have the same stride ordering.");
978 initImpl(
in.permuteStridesDescending(),
out.permuteStridesDescending(),
986 template <
class C1,
class C2>
991 vigra_precondition(
in.strideOrdering() ==
out.strideOrdering(),
992 "FFTWPlan.init(): input and output must have the same stride ordering.");
994 initImpl(
in.permuteStridesDescending(),
out.permuteStridesDescending(),
1002 template <
class C1,
class C2>
1007 vigra_precondition(
in.strideOrdering() ==
out.strideOrdering(),
1008 "FFTWPlan.init(): input and output must have the same stride ordering.");
1010 initImpl(
in.permuteStridesDescending(),
out.permuteStridesDescending(),
1021 template <
class C1,
class C2>
1025 executeImpl(
in.permuteStridesDescending(),
out.permuteStridesDescending());
1035 template <
class C1,
class C2>
1039 executeImpl(
in.permuteStridesDescending(),
out.permuteStridesDescending());
1049 template <
class C1,
class C2>
1053 executeImpl(
in.permuteStridesDescending(),
out.permuteStridesDescending());
1058 template <
class MI,
class MO>
1061 template <
class MI,
class MO>
1067 vigra_precondition(
in.shape() ==
out.shape(),
1068 "FFTWPlan.init(): input and output must have the same shape.");
1071 void checkShapes(MultiArrayView<N, Real, StridedArrayTag> ins,
1072 MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> outs)
const
1074 for(
int k=0; k<(int)N-1; ++k)
1075 vigra_precondition(ins.shape(k) == outs.shape(k),
1076 "FFTWPlan.init(): input and output must have matching shapes.");
1077 vigra_precondition(ins.shape(N-1) / 2 + 1 == outs.shape(N-1),
1078 "FFTWPlan.init(): input and output must have matching shapes.");
1081 void checkShapes(MultiArrayView<N, FFTWComplex<Real>, StridedArrayTag> ins,
1082 MultiArrayView<N, Real, StridedArrayTag> outs)
const
1084 for(
int k=0; k<(int)N-1; ++k)
1085 vigra_precondition(ins.shape(k) == outs.shape(k),
1086 "FFTWPlan.init(): input and output must have matching shapes.");
1087 vigra_precondition(outs.shape(N-1) / 2 + 1 == ins.shape(N-1),
1088 "FFTWPlan.init(): input and output must have matching shapes.");
1092template <
unsigned int N,
class Real>
1093template <
class MI,
class MO>
1095FFTWPlan<N, Real>::initImpl(MI ins, MO outs,
int SIGN,
unsigned int planner_flags)
1109 for(
unsigned int j=1;
j<N; ++
j)
1116 detail::FFTWLock<> lock;
1121 detail::fftwPlanDestroy(plan);
1131template <
unsigned int N,
class Real>
1132template <
class MI,
class MO>
1133void FFTWPlan<N, Real>::executeImpl(MI ins, MO outs)
const
1135 vigra_precondition(plan != 0,
"FFTWPlan::execute(): plan is NULL.");
1142 "FFTWPlan::execute(): shape mismatch between plan and data.");
1144 "FFTWPlan::execute(): strides mismatch between plan and input data.");
1146 "FFTWPlan::execute(): strides mismatch between plan and output data.");
1148 detail::fftwPlanExecute(plan,
ins.data(),
outs.data());
1150 typedef typename MO::value_type V;
1152 outs *= V(1.0) / Real(
outs.size());
1194template <
unsigned int N,
class Real>
1202 RArray realArray, realKernel;
1203 CArray fourierArray, fourierKernel;
1204 bool useFourierKernel;
1215 : useFourierKernel(
false)
1229 template <
class C1,
class C2,
class C3>
1234 : useFourierKernel(
false)
1250 template <
class C1,
class C2,
class C3>
1255 : useFourierKernel(
true)
1272 template <
class C1,
class C2,
class C3>
1295 template <
class C1,
class C2,
class C3>
1297 bool useFourierKernel =
false,
1300 if(useFourierKernel)
1310 template <
class C1,
class C2,
class C3>
1316 vigra_precondition(
in.shape() ==
out.shape(),
1317 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1325 template <
class C1,
class C2,
class C3>
1331 vigra_precondition(
in.shape() ==
out.shape(),
1332 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1340 template <
class C1,
class C2,
class C3>
1347 vigra_precondition(
in.shape() ==
out.shape(),
1348 "FFTWConvolvePlan::init(): input and output must have the same shape.");
1359 template <
class C1,
class KernelIterator,
class OutIterator>
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 bool realKernel = IsSameType<KernelValue, Real>::value;
1370 bool fourierKernel = IsSameType<KernelValue, Complex>::value;
1372 vigra_precondition(realKernel || fourierKernel,
1373 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
1374 vigra_precondition((IsSameType<OutValue, Real>::value),
1375 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
1384 initFourierKernelMany(
in.shape(),
1396 template <
class C1,
class KernelIterator,
class OutIterator>
1403 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1404 typedef typename KernelArray::value_type
KernelValue;
1405 typedef typename std::iterator_traits<OutIterator>::value_type
OutArray;
1406 typedef typename OutArray::value_type
OutValue;
1408 vigra_precondition((IsSameType<KernelValue, Complex>::value),
1409 "FFTWConvolvePlan::initMany(): kernels have unsuitable value_type.");
1410 vigra_precondition((IsSameType<OutValue, Complex>::value),
1411 "FFTWConvolvePlan::initMany(): outputs have unsuitable value_type.");
1422 forward_plan =
fplan;
1423 backward_plan =
bplan;
1443 void initFourierKernelMany(Shape inOut, Shape kernels,
1444 unsigned int planner_flags = FFTW_ESTIMATE)
1446 initFourierKernel(inOut, kernels, planner_flags);
1456 template <
class C1,
class C2,
class C3>
1471 template <
class C1,
class C2,
class C3>
1483 template <
class C1,
class C2,
class C3>
1496 template <
class C1,
class KernelIterator,
class OutIterator>
1508 template <
class C1,
class KernelIterator,
class OutIterator>
1513 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1514 typedef typename KernelArray::value_type
KernelValue;
1516 typedef typename std::iterator_traits<OutIterator>::value_type
OutArray;
1517 typedef typename OutArray::value_type
OutValue;
1519 bool realKernel = IsSameType<KernelValue, Real>::value;
1520 bool fourierKernel = IsSameType<KernelValue, Complex>::value;
1522 vigra_precondition(realKernel || fourierKernel,
1523 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type.");
1524 vigra_precondition((IsSameType<OutValue, Real>::value),
1525 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type.");
1532 template <
class KernelIterator,
class OutIterator>
1533 Shape checkShapes(Shape
in,
1537 template <
class KernelIterator,
class OutIterator>
1538 Shape checkShapesFourier(Shape
in,
1542 template <
class KernelIterator,
class OutIterator>
1543 Shape checkShapesComplex(Shape
in,
1547 template <
class C1,
class C2,
class C3>
1553 template <
class C1,
class KernelIterator,
class OutIterator>
1559 template <
class C1,
class KernelIterator,
class OutIterator>
1567template <
unsigned int N,
class Real>
1570 unsigned int planner_flags)
1585 forward_plan =
fplan;
1586 backward_plan =
bplan;
1591 useFourierKernel =
false;
1594template <
unsigned int N,
class Real>
1596FFTWConvolvePlan<N, Real>::initFourierKernel(Shape in, Shape kernel,
1597 unsigned int planner_flags)
1602 for(
unsigned int k=0;
k<N; ++
k)
1604 "FFTWConvolvePlan::init(): kernel too small for given input.");
1616 forward_plan =
fplan;
1617 backward_plan =
bplan;
1622 useFourierKernel =
true;
1625template <
unsigned int N,
class Real>
1627FFTWConvolvePlan<N, Real>::initComplex(Shape in, Shape kernel,
1628 unsigned int planner_flags)
1632 if(useFourierKernel)
1634 for(
unsigned int k=0;
k<N; ++
k)
1636 "FFTWConvolvePlan::init(): kernel too small for given input.");
1650 forward_plan =
fplan;
1651 backward_plan =
bplan;
1658template <
unsigned int N,
class Real>
1659template <
class C1,
class C2,
class C3>
1661FFTWConvolvePlan<N, Real>::executeImpl(MultiArrayView<N, Real, C1> in,
1662 MultiArrayView<N, Real, C2> kernel,
1663 MultiArrayView<N, Real, C3> out,
1664 bool do_correlation)
1666 vigra_precondition(!useFourierKernel,
1667 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel.");
1669 vigra_precondition(
in.shape() ==
out.shape(),
1670 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1675 right =
in.shape() + left;
1677 vigra_precondition(
paddedShape == realArray.shape(),
1678 "FFTWConvolvePlan::execute(): shape mismatch between input and plan.");
1680 detail::fftEmbedArray(
in, realArray);
1681 forward_plan.execute(realArray, fourierArray);
1683 detail::fftEmbedKernel(
kernel, realKernel);
1684 forward_plan.execute(realKernel, fourierKernel);
1688 using namespace multi_math;
1689 fourierArray *=
conj(fourierKernel);
1693 fourierArray *= fourierKernel;
1696 backward_plan.execute(fourierArray, realArray);
1698 out = realArray.subarray(left, right);
1701template <
unsigned int N,
class Real>
1702template <
class C1,
class C2,
class C3>
1705 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1706 MultiArrayView<N, Real, C3> out)
1708 vigra_precondition(useFourierKernel,
1709 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel.");
1711 vigra_precondition(
in.shape() ==
out.shape(),
1712 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1714 vigra_precondition(
kernel.shape() == fourierArray.shape(),
1715 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan.");
1720 right =
in.shape() + left;
1722 vigra_precondition(
paddedShape == realArray.shape(),
1723 "FFTWConvolvePlan::execute(): shape mismatch between input and plan.");
1725 detail::fftEmbedArray(
in, realArray);
1726 forward_plan.execute(realArray, fourierArray);
1729 moveDCToHalfspaceUpperLeft(fourierKernel);
1731 fourierArray *= fourierKernel;
1733 backward_plan.execute(fourierArray, realArray);
1735 out = realArray.subarray(left, right);
1738template <
unsigned int N,
class Real>
1739template <
class C1,
class C2,
class C3>
1742 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
1743 MultiArrayView<N, FFTWComplex<Real>, C3> out)
1745 vigra_precondition(
in.shape() ==
out.shape(),
1746 "FFTWConvolvePlan::execute(): input and output must have the same shape.");
1751 right =
in.shape() + left;
1753 if(useFourierKernel)
1755 vigra_precondition(
kernel.shape() == fourierArray.shape(),
1756 "FFTWConvolvePlan::execute(): shape mismatch between kernel and plan.");
1763 detail::fftEmbedKernel(
kernel, fourierKernel);
1764 forward_plan.execute(fourierKernel, fourierKernel);
1767 detail::fftEmbedArray(
in, fourierArray);
1768 forward_plan.execute(fourierArray, fourierArray);
1770 fourierArray *= fourierKernel;
1772 backward_plan.execute(fourierArray, fourierArray);
1774 out = fourierArray.subarray(left, right);
1777template <
unsigned int N,
class Real>
1778template <
class C1,
class KernelIterator,
class OutIterator>
1780FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
1781 KernelIterator kernels, KernelIterator kernelsEnd,
1782 OutIterator outs, VigraFalseType )
1784 vigra_precondition(!useFourierKernel,
1785 "FFTWConvolvePlan::execute(): plan was generated for Fourier kernel, got spatial kernel.");
1791 right =
in.shape() + left;
1793 vigra_precondition(
paddedShape == realArray.shape(),
1794 "FFTWConvolvePlan::executeMany(): shape mismatch between input and plan.");
1796 detail::fftEmbedArray(
in, realArray);
1797 forward_plan.execute(realArray, fourierArray);
1801 detail::fftEmbedKernel(*
kernels, realKernel);
1802 forward_plan.execute(realKernel, fourierKernel);
1804 fourierKernel *= fourierArray;
1806 backward_plan.execute(fourierKernel, realKernel);
1808 *
outs = realKernel.subarray(left, right);
1812template <
unsigned int N,
class Real>
1813template <
class C1,
class KernelIterator,
class OutIterator>
1815FFTWConvolvePlan<N, Real>::executeManyImpl(MultiArrayView<N, Real, C1> in,
1816 KernelIterator kernels, KernelIterator kernelsEnd,
1817 OutIterator outs, VigraTrueType )
1819 vigra_precondition(useFourierKernel,
1820 "FFTWConvolvePlan::execute(): plan was generated for spatial kernel, got Fourier kernel.");
1826 right =
in.shape() + left;
1828 vigra_precondition(
complexShape == fourierArray.shape(),
1829 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between kernels and plan.");
1831 vigra_precondition(
paddedShape == realArray.shape(),
1832 "FFTWConvolvePlan::executeFourierKernelMany(): shape mismatch between input and plan.");
1834 detail::fftEmbedArray(
in, realArray);
1835 forward_plan.execute(realArray, fourierArray);
1840 moveDCToHalfspaceUpperLeft(fourierKernel);
1841 fourierKernel *= fourierArray;
1843 backward_plan.execute(fourierKernel, realKernel);
1845 *
outs = realKernel.subarray(left, right);
1849template <
unsigned int N,
class Real>
1850template <
class C1,
class KernelIterator,
class OutIterator>
1853 KernelIterator kernels, KernelIterator kernelsEnd,
1856 typedef typename std::iterator_traits<KernelIterator>::value_type KernelArray;
1857 typedef typename KernelArray::value_type
KernelValue;
1858 typedef typename std::iterator_traits<OutIterator>::value_type
OutArray;
1859 typedef typename OutArray::value_type
OutValue;
1861 vigra_precondition((IsSameType<KernelValue, Complex>::value),
1862 "FFTWConvolvePlan::executeMany(): kernels have unsuitable value_type.");
1863 vigra_precondition((IsSameType<OutValue, Complex>::value),
1864 "FFTWConvolvePlan::executeMany(): outputs have unsuitable value_type.");
1869 right =
in.shape() + left;
1871 detail::fftEmbedArray(
in, fourierArray);
1872 forward_plan.execute(fourierArray, fourierArray);
1876 if(useFourierKernel)
1883 detail::fftEmbedKernel(*
kernels, fourierKernel);
1884 forward_plan.execute(fourierKernel, fourierKernel);
1887 fourierKernel *= fourierArray;
1889 backward_plan.execute(fourierKernel, fourierKernel);
1891 *
outs = fourierKernel.subarray(left, right);
1897template <
unsigned int N,
class Real>
1898template <
class KernelIterator,
class OutIterator>
1899typename FFTWConvolvePlan<N, Real>::Shape
1900FFTWConvolvePlan<N, Real>::checkShapes(Shape in,
1901 KernelIterator kernels, KernelIterator kernelsEnd,
1905 "FFTWConvolvePlan::checkShapes(): empty kernel sequence.");
1910 vigra_precondition(
in ==
outs->shape(),
1911 "FFTWConvolvePlan::checkShapes(): shape mismatch between input and (one) output.");
1915 "FFTWConvolvePlan::checkShapes(): all kernels have size 0.");
1919template <
unsigned int N,
class Real>
1920template <
class KernelIterator,
class OutIterator>
1921typename FFTWConvolvePlan<N, Real>::Shape
1922FFTWConvolvePlan<N, Real>::checkShapesFourier(Shape in,
1923 KernelIterator kernels, KernelIterator kernelsEnd,
1927 "FFTWConvolvePlan::checkShapesFourier(): empty kernel sequence.");
1932 for(
unsigned int k=0;
k<N; ++
k)
1934 "FFTWConvolvePlan::checkShapesFourier(): kernels too small for given input.");
1938 vigra_precondition(
in ==
outs->shape(),
1939 "FFTWConvolvePlan::checkShapesFourier(): shape mismatch between input and (one) output.");
1941 "FFTWConvolvePlan::checkShapesFourier(): all kernels must have the same size.");
1946template <
unsigned int N,
class Real>
1947template <
class KernelIterator,
class OutIterator>
1948typename FFTWConvolvePlan<N, Real>::Shape
1949FFTWConvolvePlan<N, Real>::checkShapesComplex(Shape in,
1950 KernelIterator kernels, KernelIterator kernelsEnd,
1954 "FFTWConvolvePlan::checkShapesComplex(): empty kernel sequence.");
1959 vigra_precondition(
in ==
outs->shape(),
1960 "FFTWConvolvePlan::checkShapesComplex(): shape mismatch between input and (one) output.");
1961 if(useFourierKernel)
1964 "FFTWConvolvePlan::checkShapesComplex(): Fourier domain kernels must have identical size.");
1972 "FFTWConvolvePlan::checkShapesComplex(): all kernels have size 0.");
1974 if(useFourierKernel)
1976 for(
unsigned int k=0;
k<N; ++
k)
1978 "FFTWConvolvePlan::checkShapesComplex(): kernels too small for given input.");
2019template <
unsigned int N,
class Real>
2047 template <
class C1,
class C2,
class C3>
2067 template <
class C1,
class C2,
class C3>
2069 bool useFourierKernel =
false,
2073 ignore_argument(useFourierKernel);
2080 template <
class C1,
class C2,
class C3>
2086 vigra_precondition(
in.shape() ==
out.shape(),
2087 "FFTWCorrelatePlan::init(): input and output must have the same shape.");
2098 template <
class C1,
class C2,
class C3>
2113template <
unsigned int N,
class Real,
class C1,
class C2>
2116 MultiArrayView<N, FFTWComplex<Real>, C2> out)
2118 FFTWPlan<N, Real>(in, out, FFTW_FORWARD).execute(in, out);
2121template <
unsigned int N,
class Real,
class C1,
class C2>
2124 MultiArrayView<N, FFTWComplex<Real>, C2> out)
2126 FFTWPlan<N, Real>(in, out, FFTW_BACKWARD).execute(in, out);
2129template <
unsigned int N,
class Real,
class C1,
class C2>
2132 MultiArrayView<N, FFTWComplex<Real>, C2> out)
2134 if(in.shape() == out.shape())
2138 FFTWPlan<N, Real>(out, out, FFTW_FORWARD).execute(out, out);
2142 FFTWPlan<N, Real>(in, out).execute(in, out);
2145 vigra_precondition(
false,
2146 "fourierTransform(): shape mismatch between input and output.");
2149template <
unsigned int N,
class Real,
class C1,
class C2>
2152 MultiArrayView<N, Real, C2> out)
2155 "fourierTransformInverse(): shape mismatch between input and output.");
2156 FFTWPlan<N, Real>(in, out).execute(in, out);
2350template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2359template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2362 MultiArrayView<N, FFTWComplex<Real>, C2> kernel,
2363 MultiArrayView<N, Real, C3> out)
2365 FFTWConvolvePlan<N, Real>(in, kernel, out).execute(in, kernel, out);
2374template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
2390template <
unsigned int N,
class Real,
class C1,
2408template <
unsigned int N,
class Real,
class C1,
2471template <
unsigned int N,
class Real,
class C1,
class C2,
class C3>
const_pointer data() const
Definition array_vector.hxx:209
Definition multi_fft.hxx:1196
FFTWConvolvePlan(Shape inOut, Shape kernel, bool useFourierKernel=false, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan from just the shape information.
Definition multi_fft.hxx:1296
void execute(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out)
Execute a plan to convolve a real array with a real kernel.
Definition multi_fft.hxx:1457
FFTWConvolvePlan(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, FFTWComplex< Real >, C3 > out, bool fourierDomainKernel, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan to convolve a complex array with a complex kernel.
Definition multi_fft.hxx:1273
void executeMany(MultiArrayView< N, Real, C1 > in, KernelIterator kernels, KernelIterator kernelsEnd, OutIterator outs)
Execute a plan to convolve a real array with a sequence of kernels.
Definition multi_fft.hxx:1509
void init(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, FFTWComplex< Real >, C3 > out, bool fourierDomainKernel, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a complex array with a complex kernel.
Definition multi_fft.hxx:1341
void init(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a real array with a complex kernel.
Definition multi_fft.hxx:1326
void execute(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, FFTWComplex< Real >, C3 > out)
Execute a plan to convolve a complex array with a complex kernel.
void execute(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, Real, C3 > out)
Execute a plan to convolve a real array with a complex kernel.
FFTWConvolvePlan(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan to convolve a real array with a real kernel.
Definition multi_fft.hxx:1230
void initMany(MultiArrayView< N, FFTWComplex< Real >, C1 > in, KernelIterator kernels, KernelIterator kernelsEnd, OutIterator outs, bool fourierDomainKernels, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a complex array with a sequence of kernels.
Definition multi_fft.hxx:1397
void executeMany(MultiArrayView< N, FFTWComplex< Real >, C1 > in, KernelIterator kernels, KernelIterator kernelsEnd, OutIterator outs)
Execute a plan to convolve a complex array with a sequence of kernels.
FFTWConvolvePlan(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan to convolve a real array with a complex kernel.
Definition multi_fft.hxx:1251
void initMany(MultiArrayView< N, Real, C1 > in, KernelIterator kernels, KernelIterator kernelsEnd, OutIterator outs, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a real array with a sequence of kernels.
Definition multi_fft.hxx:1360
FFTWConvolvePlan()
Create an empty plan.
Definition multi_fft.hxx:1214
void init(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a real array with a real kernel.
Definition multi_fft.hxx:1311
Definition multi_fft.hxx:2022
void execute(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out)
Execute a plan to correlate a real array with a real kernel.
Definition multi_fft.hxx:2099
FFTWCorrelatePlan(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan to correlate a real array with a real kernel.
Definition multi_fft.hxx:2048
void init(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, Real, C2 > kernel, MultiArrayView< N, Real, C3 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a plan to convolve a real array with a real kernel.
Definition multi_fft.hxx:2081
FFTWCorrelatePlan()
Create an empty plan.
Definition multi_fft.hxx:2032
FFTWCorrelatePlan(Shape inOut, Shape kernel, bool useFourierKernel=false, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan from just the shape information.
Definition multi_fft.hxx:2068
Definition multi_fft.hxx:853
FFTWPlan(FFTWPlan const &other)
Copy constructor.
Definition multi_fft.hxx:930
FFTWPlan(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out, int SIGN, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan for a complex-to-complex transform.
Definition multi_fft.hxx:880
~FFTWPlan()
Destructor.
Definition multi_fft.hxx:960
void init(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, Real, C2 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a complex-to-real transform.
Definition multi_fft.hxx:1003
void execute(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, Real, C2 > out) const
Execute a complex-to-real transform.
Definition multi_fft.hxx:1050
void execute(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out) const
Execute a real-to-complex transform.
Definition multi_fft.hxx:1036
FFTWPlan(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, Real, C2 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan for a complex-to-real transform.
Definition multi_fft.hxx:920
void init(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Init a real-to-complex transform.
Definition multi_fft.hxx:987
void init(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out, int SIGN, unsigned int planner_flags=FFTW_ESTIMATE)
Init a complex-to-complex transform.
Definition multi_fft.hxx:971
void execute(MultiArrayView< N, FFTWComplex< Real >, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out) const
Execute a complex-to-complex transform.
Definition multi_fft.hxx:1022
FFTWPlan & operator=(FFTWPlan const &other)
Copy assigment.
Definition multi_fft.hxx:943
FFTWPlan(MultiArrayView< N, Real, C1 > in, MultiArrayView< N, FFTWComplex< Real >, C2 > out, unsigned int planner_flags=FFTW_ESTIMATE)
Create a plan for a real-to-complex transform.
Definition multi_fft.hxx:900
FFTWPlan()
Create an empty plan.
Definition multi_fft.hxx:867
TinyVector< MultiArrayIndex, N > type
Definition multi_shape.hxx:272
Base class for, and view to, vigra::MultiArray.
Definition multi_array.hxx:705
vigra::detail::MultiIteratorChooser< StrideTag >::template Traverser< actual_dimension, T, T &, T * >::type traverser
Definition multi_array.hxx:764
void swap(MultiArray &other)
Definition multi_array.hxx:3072
Class for a single RGB value.
Definition rgbvalue.hxx:128
RGBValue()
Definition rgbvalue.hxx:209
void init(Iterator i, Iterator end)
Definition tinyvector.hxx:708
TinyVectorView< VALUETYPE, TO-FROM > subarray() const
Definition tinyvector.hxx:887
Class for fixed size vectors.
Definition tinyvector.hxx:1008
void fourierTransform(...)
Compute forward and inverse Fourier transforms.
void fourierTransformInverse(...)
Compute inverse Fourier transforms.
void convolveFFTComplex(...)
Convolve a complex-valued array by means of the Fourier transform.
void convolveFFTMany(...)
Convolve a real-valued array with a sequence of kernels by means of the Fourier transform.
TinyVector< T, N > fftwCorrespondingShapeR2C(TinyVector< T, N > shape)
Find frequency domain shape for a R2C Fourier transform.
Definition multi_fft.hxx:802
void moveDCToUpperLeft(...)
Rearrange the quadrants of a Fourier image so that the origin is in the image's upper left.
void convolveFFT(...)
Convolve an array with a kernel by means of the Fourier transform.
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition fftw3.hxx:1030
bool even(int t)
Check if an integer is even.
void correlateFFT(...)
Correlate an array with a kernel by means of the Fourier transform.
T sign(T t)
The sign function.
Definition mathutil.hxx:591
NumericTraits< V >::Promote prod(TinyVectorBase< V, SIZE, D1, D2 > const &l)
product of the vector's elements
Definition tinyvector.hxx:2097
void moveDCToCenter(...)
Rearrange the quadrants of a Fourier image so that the origin is in the image center.
std::ptrdiff_t MultiArrayIndex
Definition multi_fwd.hxx:60
FixedPoint16< IntBits3, OverflowHandling > & div(FixedPoint16< IntBits1, OverflowHandling > l, FixedPoint16< IntBits2, OverflowHandling > r, FixedPoint16< IntBits3, OverflowHandling > &result)
division with enforced result type.
Definition fixedpoint.hxx:1616
void convolveFFTComplexMany(...)
Convolve a complex-valued array with a sequence of kernels by means of the Fourier transform.
bool odd(int t)
Check if an integer is odd.
Definition metaprogramming.hxx:116