src/fft_fftw.cpp

   1 #include "fft_fftw.h"
   2
   3 #include "fft.h"
   4
   5 #ifdef OPENMP
   6   #include <omp.h>
   7 #endif
   8
   9 #if defined(ASYNC) && !defined(CUFFTW)
  10 #define FFTW_PLAN_WITH_THREADS() fftw_plan_with_nthreads(m_num_threads);
  11 #elif defined(OPENMP) && !defined(CUFFTW)
  12 #define FFTW_PLAN_WITH_THREADS() fftw_plan_with_nthreads(omp_get_max_threads());
  13 #else
  14 #define FFTW_PLAN_WITH_THREADS()
  15 #endif
  16
  17 Fftw::Fftw()
  18     : m_num_threads(1)
  19 {
  20 }
  21
  22 Fftw::Fftw(int num_threads)
  23     : m_num_threads(num_threads)
  24 {
  25 }
  26
  27 void Fftw::init(unsigned width, unsigned height, unsigned num_of_feats, unsigned num_of_scales, bool big_batch_mode)
  28 {
  29     m_width = width;
  30     m_height = height;
  31     m_num_of_feats = num_of_feats;
  32     m_num_of_scales = num_of_scales;
  33     m_big_batch_mode = big_batch_mode;
  34
  35 #if defined(ASYNC) || defined(OPENMP)
  36     fftw_init_threads();
  37 #endif //OPENMP
  38
  39 #ifndef CUFFTW
  40     std::cout << "FFT: FFTW" << std::endl;
  41 #else
  42     std::cout << "FFT: cuFFTW" << std::endl;
  43 #endif
  44     //FFT forward one scale
  45     {
  46     cv::Mat in_f = cv::Mat::zeros(m_height, m_width, CV_32FC1);
  47     ComplexMat out_f(m_height, m_width / 2 + 1, 1);
  48     plan_f = fftwf_plan_dft_r2c_2d(m_height, m_width,
  49                                    reinterpret_cast<float*>(in_f.data),
  50                                    reinterpret_cast<fftwf_complex*>(out_f.get_p_data()),
  51                                    FFTW_MEASURE);
  52     }
  53     //FFT forward all scales
  54     if(m_num_of_scales > 1 && m_big_batch_mode)
  55     {
  56     cv::Mat in_f_all = cv::Mat::zeros(m_height*m_num_of_scales, m_width, CV_32F);
  57     ComplexMat out_f_all(m_height, m_width / 2 + 1, m_num_of_scales);
  58     float *in = reinterpret_cast<float*>(in_f_all.data);
  59     fftwf_complex *out = reinterpret_cast<fftwf_complex*>(out_f_all.get_p_data());
  60     int rank = 2;
  61     int n[] = {(int)m_height, (int)m_width};
  62     int howmany = m_num_of_scales;
  63     int idist = m_height*m_width, odist = m_height*(m_width/2+1);
  64     int istride = 1, ostride = 1;
  65     int *inembed = NULL, *onembed = NULL;
  66
  67     FFTW_PLAN_WITH_THREADS();
  68     plan_f_all_scales = fftwf_plan_many_dft_r2c(rank, n, howmany,
  69                                       in, inembed, istride, idist,
  70                                       out, onembed, ostride, odist,
  71                                       FFTW_MEASURE);
  72     }
  73     //FFT forward window one scale
  74     {
  75         cv::Mat in_fw = cv::Mat::zeros(m_height * m_num_of_feats, m_width, CV_32F);
  76         ComplexMat out_fw(m_height, m_width / 2 + 1, m_num_of_feats);
  77         float *in = reinterpret_cast<float*>(in_fw.data);
  78         fftwf_complex *out = reinterpret_cast<fftwf_complex*>(out_fw.get_p_data());
  79         int rank = 2;
  80         int n[] = {(int)m_height, (int)m_width};
  81         int howmany = m_num_of_feats;
  82         int idist = m_height*m_width, odist = m_height*(m_width/2+1);
  83         int istride = 1, ostride = 1;
  84         int *inembed = NULL, *onembed = NULL;
  85
  86         FFTW_PLAN_WITH_THREADS();
  87         plan_fw = fftwf_plan_many_dft_r2c(rank, n, howmany,
  88                                           in, inembed, istride, idist,
  89                                           out, onembed, ostride, odist,
  90                                           FFTW_MEASURE);
  91     }
  92     //FFT forward window all scales all feats
  93     if(m_num_of_scales > 1 && m_big_batch_mode)
  94     {
  95         cv::Mat in_all = cv::Mat::zeros(m_height * (m_num_of_scales*m_num_of_feats), m_width, CV_32F);
  96         ComplexMat out_all(m_height, m_width / 2 + 1, m_num_of_scales*m_num_of_feats);
  97         float *in = reinterpret_cast<float*>(in_all.data);
  98         fftwf_complex *out = reinterpret_cast<fftwf_complex*>(out_all.get_p_data());
  99         int rank = 2;
 100         int n[] = {(int)m_height, (int)m_width};
 101         int howmany = m_num_of_scales*m_num_of_feats;
 102         int idist = m_height*m_width, odist = m_height*(m_width/2+1);
 103         int istride = 1, ostride = 1;
 104         int *inembed = NULL, *onembed = NULL;
 105
 106         FFTW_PLAN_WITH_THREADS();
 107         plan_fw_all_scales = fftwf_plan_many_dft_r2c(rank, n, howmany,
 108                                                      in,  inembed, istride, idist,
 109                                                      out, onembed, ostride, odist,
 110                                                      FFTW_MEASURE);
 111     }
 112     //FFT inverse one scale
 113     {
 114         ComplexMat in_i(m_height,m_width,m_num_of_feats);
 115         cv::Mat out_i = cv::Mat::zeros(m_height, m_width, CV_32FC(m_num_of_feats));
 116         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i.get_p_data());
 117         float *out = reinterpret_cast<float*>(out_i.data);
 118         int rank = 2;
 119         int n[] = {(int)m_height, (int)m_width};
 120         int howmany = m_num_of_feats;
 121         int idist = m_height*(m_width/2+1), odist = 1;
 122         int istride = 1, ostride = m_num_of_feats;
 123         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 124
 125         FFTW_PLAN_WITH_THREADS();
 126         plan_i_features = fftwf_plan_many_dft_c2r(rank, n, howmany,
 127                                                   in,  inembed, istride, idist,
 128                                                   out, onembed, ostride, odist,
 129                                                   FFTW_MEASURE);
 130     }
 131     //FFT inverse all scales
 132     if(m_num_of_scales > 1 && m_big_batch_mode)
 133     {
 134         ComplexMat in_i_all(m_height,m_width,m_num_of_feats*m_num_of_scales);
 135         cv::Mat out_i_all = cv::Mat::zeros(m_height, m_width, CV_32FC(m_num_of_feats*m_num_of_scales));
 136         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i_all.get_p_data());
 137         float *out = reinterpret_cast<float*>(out_i_all.data);
 138         int rank = 2;
 139         int n[] = {(int)m_height, (int)m_width};
 140         int howmany = m_num_of_feats*m_num_of_scales;
 141         int idist = m_height*(m_width/2+1), odist = 1;
 142         int istride = 1, ostride = m_num_of_feats*m_num_of_scales;
 143         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 144
 145         FFTW_PLAN_WITH_THREADS();
 146         plan_i_features_all_scales = fftwf_plan_many_dft_c2r(rank, n, howmany,
 147                                                   in,  inembed, istride, idist,
 148                                                   out, onembed, ostride, odist,
 149                                                   FFTW_MEASURE);
 150     }
 151     //FFT inver one channel one scale
 152     {
 153         ComplexMat in_i1(m_height,m_width,1);
 154         cv::Mat out_i1 = cv::Mat::zeros(m_height, m_width, CV_32FC1);
 155         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i1.get_p_data());
 156         float *out = reinterpret_cast<float*>(out_i1.data);
 157         int rank = 2;
 158         int n[] = {(int)m_height, (int)m_width};
 159         int howmany = 1;
 160         int idist = m_height*(m_width/2+1), odist = 1;
 161         int istride = 1, ostride = 1;
 162         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 163
 164         FFTW_PLAN_WITH_THREADS();
 165         plan_i_1ch = fftwf_plan_many_dft_c2r(rank, n, howmany,
 166                                              in,  inembed, istride, idist,
 167                                              out, onembed, ostride, odist,
 168                                              FFTW_MEASURE);
 169     }
 170     //FFT inver one channel all scales
 171     if(m_num_of_scales > 1 && m_big_batch_mode)
 172     {
 173         ComplexMat in_i1_all(m_height,m_width,m_num_of_scales);
 174         cv::Mat out_i1_all = cv::Mat::zeros(m_height, m_width, CV_32FC(m_num_of_scales));
 175         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i1_all.get_p_data());
 176         float *out = reinterpret_cast<float*>(out_i1_all.data);
 177         int rank = 2;
 178         int n[] = {(int)m_height, (int)m_width};
 179         int howmany = m_num_of_scales;
 180         int idist = m_height*(m_width/2+1), odist = 1;
 181         int istride = 1, ostride = m_num_of_scales;
 182         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 183
 184         FFTW_PLAN_WITH_THREADS();
 185         plan_i_1ch_all_scales = fftwf_plan_many_dft_c2r(rank, n, howmany,
 186                                              in,  inembed, istride, idist,
 187                                              out, onembed, ostride, odist,
 188                                              FFTW_MEASURE);
 189     }
 190 }
 191
 192 void Fftw::set_window(const cv::Mat &window)
 193 {
 194     m_window = window;
 195 }
 196
 197 ComplexMat Fftw::forward(const cv::Mat &input)
 198 {
 199     ComplexMat complex_result;
 200     if(m_big_batch_mode && input.rows == (int)(m_height*m_num_of_scales)){
 201         complex_result.create(m_height, m_width / 2 + 1, m_num_of_scales);
 202         fftwf_execute_dft_r2c(plan_f_all_scales, reinterpret_cast<float*>(input.data),
 203                                 reinterpret_cast<fftwf_complex*>(complex_result.get_p_data()));
 204     } else {
 205         complex_result.create(m_height, m_width / 2 + 1, 1);
 206         fftwf_execute_dft_r2c(plan_f, reinterpret_cast<float*>(input.data),
 207                                 reinterpret_cast<fftwf_complex*>(complex_result.get_p_data()));
 208     }
 209     return complex_result;
 210 }
 211
 212 ComplexMat Fftw::forward_window(const std::vector<cv::Mat> &input)
 213 {
 214     int n_channels = input.size();
 215     cv::Mat in_all(m_height * n_channels, m_width, CV_32F);
 216     for (int i = 0; i < n_channels; ++i) {
 217         cv::Mat in_roi(in_all, cv::Rect(0, i*m_height, m_width, m_height));
 218         in_roi = input[i].mul(m_window);
 219     }
 220     ComplexMat result;
 221     if(n_channels > (int) m_num_of_feats)
 222         result.create(m_height, m_width/2 + 1, n_channels,m_num_of_scales);
 223     else
 224         result.create(m_height, m_width/2 + 1, n_channels);
 225
 226     float *in = reinterpret_cast<float*>(in_all.data);
 227     fftwf_complex *out = reinterpret_cast<fftwf_complex*>(result.get_p_data());
 228
 229     if (n_channels <= (int) m_num_of_feats)
 230         fftwf_execute_dft_r2c(plan_fw, in, out);
 231     else
 232         fftwf_execute_dft_r2c(plan_fw_all_scales, in, out);
 233
 234     return result;
 235 }
 236
 237 cv::Mat Fftw::inverse(const ComplexMat &inputf)
 238 {
 239     int n_channels = inputf.n_channels;
 240     cv::Mat real_result(m_height, m_width, CV_32FC(n_channels));
 241     fftwf_complex *in = reinterpret_cast<fftwf_complex*>(inputf.get_p_data());
 242     float *out = reinterpret_cast<float*>(real_result.data);
 243
 244     if(n_channels == 1)
 245         fftwf_execute_dft_c2r(plan_i_1ch, in, out);
 246     else if(m_big_batch_mode && n_channels == (int) m_num_of_scales)
 247         fftwf_execute_dft_c2r(plan_i_1ch_all_scales, in, out);
 248     else if(m_big_batch_mode && n_channels == (int) m_num_of_feats * (int) m_num_of_scales)
 249         fftwf_execute_dft_c2r(plan_i_features_all_scales, in, out);
 250     else
 251         fftwf_execute_dft_c2r(plan_i_features, in, out);
 252
 253     return real_result/(m_width*m_height);
 254 }
 255
 256 Fftw::~Fftw()
 257 {
 258   fftwf_destroy_plan(plan_f);
 259   fftwf_destroy_plan(plan_f_all_scales);
 260   fftwf_destroy_plan(plan_fw);
 261   fftwf_destroy_plan(plan_fw_all_scales);
 262   fftwf_destroy_plan(plan_i_features);
 263   fftwf_destroy_plan(plan_i_features_all_scales);
 264   fftwf_destroy_plan(plan_i_1ch);
 265   fftwf_destroy_plan(plan_i_1ch_all_scales);
 266 }