SourceXtractorPlusPlus  0.13
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FFT.cpp
Go to the documentation of this file.
1 
17 /*
18  * @file FFT.cpp
19  * @date 11/09/18
20  * @author aalvarez
21  */
22 
23 #include <fftw3.h>
24 #include <map>
25 #include <boost/thread/shared_mutex.hpp>
26 #include "SEFramework/FFT/FFT.h"
27 
28 namespace SourceXtractor {
29 
34 boost::mutex fftw_global_plan_mutex {};
35 
41 
47 
48 
49 int fftRoundDimension(int size) {
50  // Precomputed lookup table for the optimal dimension
51  static int optimal_dims[] = {
52  1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
53  11, 12, 13, 14, 15, 16, 18, 18, 20, 20, 21,
54  22, 24, 24, 25, 26, 27, 28, 30, 30, 32, 32,
55  33, 35, 35, 36, 39, 39, 39, 40, 42, 42, 44,
56  44, 45, 48, 48, 48, 49, 50, 52, 52, 54, 54,
57  55, 56, 60, 60, 60, 60, 63, 63, 63, 64, 65,
58  66, 70, 70, 70, 70, 72, 72, 75, 75, 75, 77,
59  77, 78, 80, 80, 81, 84, 84, 84, 88, 88, 88,
60  88, 90, 90, 91, 96, 96, 96, 96, 96, 98, 98,
61  99, 100, 104, 104, 104, 104, 105, 108, 108, 108, 110,
62  110, 112, 112, 117, 117, 117, 117, 117, 120, 120, 120,
63  125, 125, 125, 125, 125, 126, 128, 128, 130, 130, 132,
64  132, 135, 135, 135, 140, 140, 140, 140, 140, 144, 144,
65  144, 144, 147, 147, 147, 150, 150, 150, 154, 154, 154,
66  154, 156, 156, 160, 160, 160, 160, 162, 162, 165, 165,
67  165, 168, 168, 168, 175, 175, 175, 175, 175, 175, 175,
68  176, 180, 180, 180, 180, 182, 182, 189, 189, 189, 189,
69  189, 189, 189, 192, 192, 192, 195, 195, 195, 196, 198,
70  198, 200, 200, 208, 208, 208, 208, 208, 208, 208, 208,
71  210, 210, 216, 216, 216, 216, 216, 216, 220, 220, 220,
72  220, 224, 224, 224, 224, 225, 231, 231, 231, 231, 231,
73  231, 234, 234, 234, 240, 240, 240, 240, 240, 240, 243,
74  243, 243, 245, 245, 250, 250, 250, 250, 250, 252, 252,
75  256, 256, 256, 256, 260, 260, 260, 260, 264, 264, 264,
76  264, 270, 270, 270, 270, 270, 270, 273, 273, 273, 275,
77  275, 280, 280, 280, 280, 280, 288, 288, 288, 288, 288,
78  288, 288, 288, 294, 294, 294, 294, 294, 294, 297, 297,
79  297, 300, 300, 300, 308, 308, 308, 308, 308, 308, 308,
80  308, 312, 312, 312, 312, 315, 315, 315, 320, 320, 320,
81  320, 320, 324, 324, 324, 324, 325, 330, 330, 330, 330,
82  330, 336, 336, 336, 336, 336, 336, 343, 343, 343, 343,
83  343, 343, 343, 350, 350, 350, 350, 350, 350, 350, 351,
84  352, 360, 360, 360, 360, 360, 360, 360, 360, 364, 364,
85  364, 364, 375, 375, 375, 375, 375, 375, 375, 375, 375,
86  375, 375, 378, 378, 378, 384, 384, 384, 384, 384, 384,
87  385, 390, 390, 390, 390, 390, 392, 392, 396, 396, 396,
88  396, 400, 400, 400, 400, 405, 405, 405, 405, 405, 416,
89  416, 416, 416, 416, 416, 416, 416, 416, 416, 416, 420,
90  420, 420, 420, 432, 432, 432, 432, 432, 432, 432, 432,
91  432, 432, 432, 432, 440, 440, 440, 440, 440, 440, 440,
92  440, 441, 448, 448, 448, 448, 448, 448, 448, 450, 450,
93  455, 455, 455, 455, 455, 462, 462, 462, 462, 462, 462,
94  462, 468, 468, 468, 468, 468, 468, 480, 480, 480, 480,
95  480, 480, 480, 480, 480, 480, 480, 480, 486, 486, 486,
96  486, 486, 486, 490, 490, 490, 490, 495, 495, 495, 495,
97  495, 500, 500, 500, 500, 500, 504, 504, 504, 504, 512,
98  512, 512, 512, 512, 512, 512, 512, 520, 520, 520, 520,
99  520, 520, 520, 520, 525, 525, 525, 525, 525, 528, 528,
100  528, 539, 539, 539, 539, 539, 539, 539, 539, 539, 539,
101  539, 540, 546, 546, 546, 546, 546, 546, 550, 550, 550,
102  550, 560, 560, 560, 560, 560, 560, 560, 560, 560, 560,
103  567, 567, 567, 567, 567, 567, 567, 576, 576, 576, 576,
104  576, 576, 576, 576, 576, 585, 585, 585, 585, 585, 585,
105  585, 585, 585, 588, 588, 588, 594, 594, 594, 594, 594,
106  594, 600, 600, 600, 600, 600, 600, 616, 616, 616, 616,
107  616, 616, 616, 616, 616, 616, 616, 616, 616, 616, 616,
108  616, 624, 624, 624, 624, 624, 624, 624, 624, 625, 630,
109  630, 630, 630, 630, 637, 637, 637, 637, 637, 637, 637,
110  640, 640, 640, 648, 648, 648, 648, 648, 648, 648, 648,
111  650, 650, 660, 660, 660, 660, 660, 660, 660, 660, 660,
112  660, 672, 672, 672, 672, 672, 672, 672, 672, 672, 672,
113  672, 672, 675, 675, 675, 686, 686, 686, 686, 686, 686,
114  686, 686, 686, 686, 686, 693, 693, 693, 693, 693, 693,
115  693, 700, 700, 700, 700, 700, 700, 700, 702, 702, 704,
116  704, 720, 720, 720, 720, 720, 720, 720, 720, 720, 720,
117  720, 720, 720, 720, 720, 720, 728, 728, 728, 728, 728,
118  728, 728, 728, 729, 735, 735, 735, 735, 735, 735, 750,
119  750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750,
120  750, 750, 750, 756, 756, 756, 756, 756, 756, 768, 768,
121  768, 768, 768, 768, 768, 768, 768, 768, 768, 768, 770,
122  770, 780, 780, 780, 780, 780, 780, 780, 780, 780, 780,
123  784, 784, 784, 784, 792, 792, 792, 792, 792, 792, 792,
124  792, 800, 800, 800, 800, 800, 800, 800, 800, 810, 810,
125  810, 810, 810, 810, 810, 810, 810, 810, 819, 819, 819,
126  819, 819, 819, 819, 819, 819, 825, 825, 825, 825, 825,
127  825, 832, 832, 832, 832, 832, 832, 832, 840, 840, 840,
128  840, 840, 840, 840, 840, 864, 864, 864, 864, 864, 864,
129  864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864,
130  864, 864, 864, 864, 864, 864, 864, 875, 875, 875, 875,
131  875, 875, 875, 875, 875, 875, 875, 880, 880, 880, 880,
132  880, 882, 882, 891, 891, 891, 891, 891, 891, 891, 891,
133  891, 896, 896, 896, 896, 896, 900, 900, 900, 900, 910,
134  910, 910, 910, 910, 910, 910, 910, 910, 910, 924, 924,
135  924, 924, 924, 924, 924, 924, 924, 924, 924, 924, 924,
136  924, 936, 936, 936, 936, 936, 936, 936, 936, 936, 936,
137  936, 936, 945, 945, 945, 945, 945, 945, 945, 945, 945,
138  960, 960, 960, 960, 960, 960, 960, 960, 960, 960, 960,
139  960, 960, 960, 960, 972, 972, 972, 972, 972, 972, 972,
140  972, 972, 972, 972, 972, 975, 975, 975, 980, 980, 980,
141  980, 980, 990, 990, 990, 990, 990, 990, 990, 990, 990,
142  990, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
143  1008, 1008, 1008, 1008, 1008, 1008, 1008, 1008, 1024, 1024, 1024,
144  };
145  static int n_optimal_dims = sizeof(optimal_dims) / sizeof(int);
146 
147  if (size < n_optimal_dims) {
148  return optimal_dims[size];
149  }
150  return (size / 512 + (size % 512 != 0)) * 512;
151 }
152 
153 
154 template <typename T>
156  // Make sure the buffers are big enough
157  assert(in.size() >= static_cast<size_t>(width * height * howmany));
158  assert(out.size() >= static_cast<size_t>(width * height * howmany));
159 
160  // Cache plan, as they can be reused
161  static boost::shared_mutex mutex;
163 
164  boost::upgrade_lock<boost::shared_mutex> read_lock{mutex};
165 
166  auto pi = plan_cache.find(std::make_tuple(howmany, width, height));
167  if (pi != plan_cache.end()) {
168  return pi->second;
169  }
170 
171  // No available plan yet, so get one from FFTW
172  boost::upgrade_to_unique_lock<boost::shared_mutex> write_lock{read_lock};
173  boost::lock_guard<boost::mutex> lock_planner{fftw_global_plan_mutex};
174 
175  int dims[] = {height, width};
176  int total_size = height * width;
177  int rank = 2; // 2D only
178  int istride = 1, ostride = 1;
179  int idist = total_size;
180  int odist = total_size;
181 
182  pi = plan_cache.emplace(
183  std::make_tuple(howmany, width, height),
184  plan_ptr_t{
185  fftw_traits::func_plan_fwd(
186  rank, dims, howmany,
187  in.data(), nullptr, istride, idist,
188  reinterpret_cast<typename fftw_traits::complex_t*>(out.data()), nullptr, ostride, odist,
189  FFTW_ESTIMATE | FFTW_DESTROY_INPUT
190  ),
191  fftw_traits::func_destroy_plan
192  }
193  ).first;
194 
195  return pi->second;
196 }
197 
198 template <typename T>
200  // Make sure the buffers are big enough
201  assert(in.size() >= static_cast<size_t>(width * height * howmany));
202  assert(out.size() >= static_cast<size_t>(width * height * howmany));
203 
204  // Cache plan, as they can be reused
205  static boost::shared_mutex mutex;
207 
208  boost::upgrade_lock<boost::shared_mutex> read_lock{mutex};
209 
210  auto pi = plan_cache.find(std::make_tuple(howmany, width, height));
211  if (pi != plan_cache.end()) {
212  return pi->second;
213  }
214 
215  // No available plan yet, so get one from FFTW
216  boost::upgrade_to_unique_lock<boost::shared_mutex> write_lock{read_lock};
217  boost::lock_guard<boost::mutex> lock_planner{fftw_global_plan_mutex};
218 
219  int dims[] = {height, width};
220  int total_size = height * width;
221  int rank = 2;
222  int istride = 1, ostride = 1;
223  int idist = total_size;
224  int odist = total_size;
225 
226  pi = plan_cache.emplace(
227  std::make_tuple(howmany, width, height),
228  plan_ptr_t{
229  fftw_traits::func_plan_inv(
230  rank, dims, howmany,
231  reinterpret_cast<typename fftw_traits::complex_t*>(in.data()), nullptr, istride, idist,
232  out.data(), nullptr, ostride, odist,
233  FFTW_ESTIMATE | FFTW_DESTROY_INPUT
234  ),
235  fftw_traits::func_destroy_plan
236  }
237  ).first;
238 
239  return pi->second;
240 }
241 
242 
243 template<typename T>
245  fftw_traits::func_execute_fwd(plan.get(), in.data(), reinterpret_cast<typename fftw_traits::complex_t *>(out.data()));
246 }
247 
248 
249 template<typename T>
251  fftw_traits::func_execute_inv(plan.get(), reinterpret_cast<typename fftw_traits::complex_t *>(in.data()), out.data());
252 }
253 
254 
255 template class FFT<float>;
256 template class FFT<double>;
257 
258 } // end SourceXtractor
constexpr double pi
static plan_ptr_t createForwardPlan(int howmany, int width, int height, std::vector< T > &in, std::vector< complex_t > &out)
Definition: FFT.cpp:155
static func_execute_inv_t * func_execute_inv
Definition: FFT.h:57
static func_plan_inv_t * func_plan_inv
Definition: FFT.h:74
static func_plan_fwd_t * func_plan_fwd
Definition: FFT.h:53
static plan_ptr_t createInversePlan(int howmany, int width, int height, std::vector< complex_t > &in, std::vector< T > &out)
Definition: FFT.cpp:199
T make_tuple(T...args)
T end(T...args)
STL class.
decltype(fftwf_execute_dft_r2c) typedef func_execute_fwd_t
Definition: FFT.h:50
Wraps the FFTW API with a more C++ like one.
Definition: FFT.h:87
static func_execute_fwd_t * func_execute_fwd
Definition: FFT.h:56
static void executeInverse(plan_ptr_t &plan, std::vector< complex_t > &in, std::vector< T > &out)
Definition: FFT.cpp:250
decltype(fftw_plan_many_dft_c2r) typedef func_plan_inv_t
Definition: FFT.h:68
T data(T...args)
decltype(fftw_plan_many_dft_r2c) typedef func_plan_fwd_t
Definition: FFT.h:67
decltype(fftwf_execute_dft_c2r) typedef func_execute_inv_t
Definition: FFT.h:51
static func_plan_fwd_t * func_plan_fwd
Definition: FFT.h:73
static func_execute_inv_t * func_execute_inv
Definition: FFT.h:77
decltype(fftw_execute_dft_r2c) typedef func_execute_fwd_t
Definition: FFT.h:70
decltype(fftw_destroy_plan) typedef func_destroy_plan_t
Definition: FFT.h:69
decltype(fftwf_plan_many_dft_c2r) typedef func_plan_inv_t
Definition: FFT.h:48
T get(T...args)
decltype(fftw_execute_dft_c2r) typedef func_execute_inv_t
Definition: FFT.h:71
T find(T...args)
static func_execute_fwd_t * func_execute_fwd
Definition: FFT.h:76
STL class.
int fftRoundDimension(int size)
Definition: FFT.cpp:49
decltype(fftwf_destroy_plan) typedef func_destroy_plan_t
Definition: FFT.h:49
decltype(fftwf_plan_many_dft_r2c) typedef func_plan_fwd_t
Definition: FFT.h:47
boost::mutex fftw_global_plan_mutex
Definition: FFT.cpp:34
T emplace(T...args)
static func_destroy_plan_t * func_destroy_plan
Definition: FFT.h:75
static func_destroy_plan_t * func_destroy_plan
Definition: FFT.h:55
static void executeForward(plan_ptr_t &plan, std::vector< T > &in, std::vector< complex_t > &out)
Definition: FFT.cpp:244
static func_plan_inv_t * func_plan_inv
Definition: FFT.h:54