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