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(2)
  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     //FFT forward all scales
  52     if (m_num_of_scales > 1 && m_big_batch_mode) {
  53         cv::Mat in_f_all = cv::Mat::zeros(m_height*m_num_of_scales, m_width, CV_32F);
  54         ComplexMat out_f_all(m_height, m_width / 2 + 1, m_num_of_scales);
  55         float *in = reinterpret_cast<float*>(in_f_all.data);
  56         fftwf_complex *out = reinterpret_cast<fftwf_complex*>(out_f_all.get_p_data());
  57         int rank = 2;
  58         int n[] = {(int)m_height, (int)m_width};
  59         int howmany = m_num_of_scales;
  60         int idist = m_height*m_width, odist = m_height*(m_width/2+1);
  61         int istride = 1, ostride = 1;
  62         int *inembed = NULL, *onembed = NULL;
  63
  64         FFTW_PLAN_WITH_THREADS();
  65         plan_f_all_scales = fftwf_plan_many_dft_r2c(rank, n, howmany,
  66                                                     in, inembed, istride, idist,
  67                                                     out, onembed, ostride, odist,
  68                                                     FFTW_PATIENT);
  69     }
  70     //FFT forward window one scale
  71     {
  72         cv::Mat in_fw = cv::Mat::zeros(m_height * m_num_of_feats, m_width, CV_32F);
  73         ComplexMat out_fw(m_height, m_width / 2 + 1, m_num_of_feats);
  74         float *in = reinterpret_cast<float*>(in_fw.data);
  75         fftwf_complex *out = reinterpret_cast<fftwf_complex*>(out_fw.get_p_data());
  76         int rank = 2;
  77         int n[] = {(int)m_height, (int)m_width};
  78         int howmany = m_num_of_feats;
  79         int idist = m_height*m_width, odist = m_height*(m_width/2+1);
  80         int istride = 1, ostride = 1;
  81         int *inembed = NULL, *onembed = NULL;
  82
  83         FFTW_PLAN_WITH_THREADS();
  84         plan_fw = fftwf_plan_many_dft_r2c(rank, n, howmany,
  85                                           in, inembed, istride, idist,
  86                                           out, onembed, ostride, odist,
  87                                           FFTW_PATIENT);
  88     }
  89     //FFT forward window all scales all feats
  90     if (m_num_of_scales > 1 && m_big_batch_mode) {
  91         cv::Mat in_all = cv::Mat::zeros(m_height * (m_num_of_scales*m_num_of_feats), m_width, CV_32F);
  92         ComplexMat out_all(m_height, m_width / 2 + 1, m_num_of_scales*m_num_of_feats);
  93         float *in = reinterpret_cast<float*>(in_all.data);
  94         fftwf_complex *out = reinterpret_cast<fftwf_complex*>(out_all.get_p_data());
  95         int rank = 2;
  96         int n[] = {(int)m_height, (int)m_width};
  97         int howmany = m_num_of_scales*m_num_of_feats;
  98         int idist = m_height*m_width, odist = m_height*(m_width/2+1);
  99         int istride = 1, ostride = 1;
 100         int *inembed = NULL, *onembed = NULL;
 101
 102         FFTW_PLAN_WITH_THREADS();
 103         plan_fw_all_scales = fftwf_plan_many_dft_r2c(rank, n, howmany,
 104                                                      in,  inembed, istride, idist,
 105                                                      out, onembed, ostride, odist,
 106                                                      FFTW_PATIENT);
 107     }
 108     //FFT inverse one scale
 109     {
 110         ComplexMat in_i(m_height,m_width,m_num_of_feats);
 111         cv::Mat out_i = cv::Mat::zeros(m_height, m_width, CV_32FC(m_num_of_feats));
 112         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i.get_p_data());
 113         float *out = reinterpret_cast<float*>(out_i.data);
 114         int rank = 2;
 115         int n[] = {(int)m_height, (int)m_width};
 116         int howmany = m_num_of_feats;
 117         int idist = m_height*(m_width/2+1), odist = 1;
 118         int istride = 1, ostride = m_num_of_feats;
 119         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 120
 121         FFTW_PLAN_WITH_THREADS();
 122         plan_i_features = fftwf_plan_many_dft_c2r(rank, n, howmany,
 123                                                   in,  inembed, istride, idist,
 124                                                   out, onembed, ostride, odist,
 125                                                   FFTW_PATIENT);
 126     }
 127     //FFT inverse all scales
 128     if (m_num_of_scales > 1 && m_big_batch_mode) {
 129         ComplexMat in_i_all(m_height,m_width,m_num_of_feats*m_num_of_scales);
 130         cv::Mat out_i_all = cv::Mat::zeros(m_height, m_width, CV_32FC(m_num_of_feats*m_num_of_scales));
 131         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i_all.get_p_data());
 132         float *out = reinterpret_cast<float*>(out_i_all.data);
 133         int rank = 2;
 134         int n[] = {(int)m_height, (int)m_width};
 135         int howmany = m_num_of_feats*m_num_of_scales;
 136         int idist = m_height*(m_width/2+1), odist = 1;
 137         int istride = 1, ostride = m_num_of_feats*m_num_of_scales;
 138         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 139
 140         FFTW_PLAN_WITH_THREADS();
 141         plan_i_features_all_scales = fftwf_plan_many_dft_c2r(rank, n, howmany,
 142                                                              in,  inembed, istride, idist,
 143                                                              out, onembed, ostride, odist,
 144                                                              FFTW_PATIENT);
 145     }
 146     //FFT inver one channel one scale
 147     if(!m_big_batch_mode) {
 148         ComplexMat in_i1(m_height,m_width,1);
 149         cv::Mat out_i1 = cv::Mat::zeros(m_height, m_width, CV_32FC1);
 150         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i1.get_p_data());
 151         float *out = reinterpret_cast<float*>(out_i1.data);
 152         int rank = 2;
 153         int n[] = {(int)m_height, (int)m_width};
 154         int howmany = 1;
 155         int idist = m_height*(m_width/2+1), odist = 1;
 156         int istride = 1, ostride = 1;
 157         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 158
 159         FFTW_PLAN_WITH_THREADS();
 160         plan_i_1ch = fftwf_plan_many_dft_c2r(rank, n, howmany,
 161                                              in,  inembed, istride, idist,
 162                                              out, onembed, ostride, odist,
 163                                              FFTW_PATIENT);
 164     }
 165     //FFT inver one channel all scales
 166     if (m_num_of_scales > 1 && m_big_batch_mode) {
 167         ComplexMat in_i1_all(m_height,m_width,m_num_of_scales);
 168         cv::Mat out_i1_all = cv::Mat::zeros(m_height, m_width, CV_32FC(m_num_of_scales));
 169         fftwf_complex *in = reinterpret_cast<fftwf_complex*>(in_i1_all.get_p_data());
 170         float *out = reinterpret_cast<float*>(out_i1_all.data);
 171         int rank = 2;
 172         int n[] = {(int)m_height, (int)m_width};
 173         int howmany = m_num_of_scales;
 174         int idist = m_height*(m_width/2+1), odist = 1;
 175         int istride = 1, ostride = m_num_of_scales;
 176         int inembed[] = {(int)m_height, (int)m_width/2+1}, *onembed = n;
 177
 178         FFTW_PLAN_WITH_THREADS();
 179         plan_i_1ch_all_scales = fftwf_plan_many_dft_c2r(rank, n, howmany,
 180                                                         in,  inembed, istride, idist,
 181                                                         out, onembed, ostride, odist,
 182                                                         FFTW_PATIENT);
 183     }
 184 }
 185
 186 void Fftw::set_window(const cv::Mat &window)
 187 {
 188     m_window = window;
 189 }
 190
 191 ComplexMat Fftw::forward(const cv::Mat & input)
 192 {
 193     ComplexMat complex_result;
 194     if(m_big_batch_mode && input.rows == (int)(m_height*m_num_of_scales)){
 195         complex_result.create(m_height, m_width / 2 + 1, m_num_of_scales);
 196         fftwf_execute_dft_r2c(plan_f_all_scales, reinterpret_cast<float*>(input.data),
 197                               reinterpret_cast<fftwf_complex*>(complex_result.get_p_data()));
 198     } else {
 199         complex_result.create(m_height, m_width / 2 + 1, 1);
 200         fftwf_execute_dft_r2c(plan_f, reinterpret_cast<float*>(input.data),
 201                               reinterpret_cast<fftwf_complex*>(complex_result.get_p_data()));
 202     }
 203     return complex_result;
 204 }
 205
 206 void Fftw::forward(Scale_vars & vars)
 207 {
 208     ComplexMat *complex_result = vars.flag & Track_flags::AUTO_CORRELATION ? & vars.kf : & vars.kzf;
 209     if(m_big_batch_mode && vars.in_all.rows == (int)(m_height*m_num_of_scales)){
 210         fftwf_execute_dft_r2c(plan_f_all_scales, reinterpret_cast<float*>(vars.in_all.data),
 211                               reinterpret_cast<fftwf_complex*>(complex_result->get_p_data()));
 212     } else {
 213         fftwf_execute_dft_r2c(plan_f, reinterpret_cast<float*>(vars.in_all.data),
 214                               reinterpret_cast<fftwf_complex*>(complex_result->get_p_data()));
 215     }
 216     return;
 217 }
 218
 219 ComplexMat Fftw::forward_raw(float *input, bool all_scales)
 220 {
 221     ComplexMat dummy;
 222     return dummy;
 223 }
 224
 225 ComplexMat Fftw::forward_window(const std::vector<cv::Mat> & input)
 226 {
 227     int n_channels = input.size();
 228     cv::Mat in_all(m_height * n_channels, m_width, CV_32F);
 229     for (int i = 0; i < n_channels; ++i) {
 230         cv::Mat in_roi(in_all, cv::Rect(0, i*m_height, m_width, m_height));
 231         in_roi = input[i].mul(m_window);
 232     }
 233     ComplexMat result;
 234     if(n_channels > (int) m_num_of_feats)
 235         result.create(m_height, m_width/2 + 1, n_channels,m_num_of_scales);
 236     else
 237         result.create(m_height, m_width/2 + 1, n_channels);
 238
 239     float *in = reinterpret_cast<float*>(in_all.data);
 240     fftwf_complex *out = reinterpret_cast<fftwf_complex*>(result.get_p_data());
 241
 242     if (n_channels <= (int) m_num_of_feats)
 243         fftwf_execute_dft_r2c(plan_fw, in, out);
 244     else
 245         fftwf_execute_dft_r2c(plan_fw_all_scales, in, out);
 246
 247     return result;
 248 }
 249
 250 void Fftw::forward_window(Scale_vars & vars)
 251 {
 252     int n_channels = vars.patch_feats.size();
 253     cv::Mat in_all(m_height * n_channels, m_width, CV_32F);
 254     for (int i = 0; i < n_channels; ++i) {
 255         cv::Mat in_roi(in_all, cv::Rect(0, i*m_height, m_width, m_height));
 256         in_roi = vars.patch_feats[i].mul(m_window);
 257     }
 258     ComplexMat *result = vars.flag & Track_flags::TRACKER_UPDATE ? & vars.xf : & vars.zf;
 259
 260     float *in = reinterpret_cast<float*>(in_all.data);
 261     fftwf_complex *out = reinterpret_cast<fftwf_complex*>(result->get_p_data());
 262
 263     if (n_channels <= (int) m_num_of_feats)
 264         fftwf_execute_dft_r2c(plan_fw, in, out);
 265     else
 266         fftwf_execute_dft_r2c(plan_fw_all_scales, in, out);
 267     return;
 268 }
 269
 270 cv::Mat Fftw::inverse(const ComplexMat &input)
 271 {
 272     int n_channels = input.n_channels;
 273     cv::Mat real_result(m_height, m_width, CV_32FC(n_channels));
 274     fftwf_complex *in = reinterpret_cast<fftwf_complex*>(input.get_p_data());
 275     float *out = reinterpret_cast<float*>(real_result.data);
 276
 277     if(n_channels == 1)
 278         fftwf_execute_dft_c2r(plan_i_1ch, in, out);
 279     else if(m_big_batch_mode && n_channels == (int) m_num_of_scales)
 280         fftwf_execute_dft_c2r(plan_i_1ch_all_scales, in, out);
 281     else if(m_big_batch_mode && n_channels == (int) m_num_of_feats * (int) m_num_of_scales)
 282         fftwf_execute_dft_c2r(plan_i_features_all_scales, in, out);
 283     else
 284         fftwf_execute_dft_c2r(plan_i_features, in, out);
 285
 286     return real_result/(m_width*m_height);
 287 }
 288
 289 void Fftw::inverse(Scale_vars & vars)
 290 {
 291     ComplexMat *input = vars.flag & Track_flags::RESPONSE ? & vars.kzf : &  vars.xyf;
 292     cv::Mat *real_result = vars.flag & Track_flags::RESPONSE ? & vars.response : & vars.ifft2_res;
 293
 294     int n_channels = input->n_channels;
 295     fftwf_complex *in = reinterpret_cast<fftwf_complex*>(input->get_p_data());
 296     float *out = reinterpret_cast<float*>(real_result->data);
 297
 298     if(n_channels == 1)
 299         fftwf_execute_dft_c2r(plan_i_1ch, in, out);
 300     else if(m_big_batch_mode && n_channels == (int) m_num_of_scales)
 301         fftwf_execute_dft_c2r(plan_i_1ch_all_scales, in, out);
 302     else if(m_big_batch_mode && n_channels == (int) m_num_of_feats * (int) m_num_of_scales)
 303         fftwf_execute_dft_c2r(plan_i_features_all_scales, in, out);
 304     else
 305         fftwf_execute_dft_c2r(plan_i_features, in, out);
 306
 307     *real_result = *real_result/(m_width*m_height);
 308     return;
 309 }
 310
 311 float* Fftw::inverse_raw(const ComplexMat &input)
 312 {
 313     return nullptr;
 314 }
 315
 316 Fftw::~Fftw()
 317 {
 318     fftwf_destroy_plan(plan_f);
 319     fftwf_destroy_plan(plan_fw);
 320     fftwf_destroy_plan(plan_i_features);
 321     fftwf_destroy_plan(plan_i_1ch);
 322
 323     if (m_big_batch_mode) {
 324         fftwf_destroy_plan(plan_f_all_scales);
 325         fftwf_destroy_plan(plan_i_features_all_scales);
 326         fftwf_destroy_plan(plan_fw_all_scales);
 327         fftwf_destroy_plan(plan_i_1ch_all_scales);
 328     }
 329 }