Inverse Flow - Forward Warping or Bilinear "Splatting"

edit

I am interested in the inverse or backward flow given the forward flow. The following code does this forward warping or bilinear splatting, but it is annoyingly slow (~8ms @VGA on my i7-7820HK). It seems likely to me that this could/should be closer to 1-2ms. Any insights into speeding this up?

inline bool isOnImage(const cv::Point& pt, const cv::Size& size)
{
  return pt.x >= 0 && pt.x < size.width && pt.y >= 0 && pt.y < size.height;    
}

cv::Mat img_proc::inverseFlow(const cv::Mat& flow)
{ 
  cv::Mat inverse_flow = cv::Mat::zeros(flow.size(), CV_32FC2);
  cv::Mat weights = cv::Mat::zeros(flow.size(), CV_32FC2);

  const int rows = flow.rows;
  const int cols = flow.cols;

  for(int i = 0; i < rows; ++i)
  {
    auto flow_ptr = flow.ptr<cv::Vec2f>(i);

    for(int j = 0; j < cols; ++j)
    {
      const float du = flow_ptr[j][0];
      const float dv = flow_ptr[j][1];

      const int u = j + std::round(du);
      const int v = i + std::round(dv);                                         

      if(!isOnImage({u,v}, flow.size()))                                        
      {                
        continue;                       
      }

      const int du_floor = (int) std::floor(du);
      const int du_ceil = (int) std::ceil(du);   

      const int dv_floor = (int) std::floor(dv);
      const int dv_ceil = (int) std::ceil(dv);

      const int u_min = std::min(cols-1, std::max(0, j + du_floor));
      const int u_max = std::min(cols-1, std::max(0, j + du_ceil));

      const int v_min = std::min(rows-1, std::max(0, i + dv_floor));
      const int v_max = std::min(rows-1, std::max(0, i + dv_ceil));

      const float uf = j + du;
      const float vf = i + dv;

      const float w0 = (u_max - uf) * (v_max - vf); // TL
      const float w1 = (uf - u_min) * (v_max - vf); // TR
      const float w2 = (uf - u_min) * (vf - v_min); // BR
      const float w3 = (u_max - uf) * (vf - v_min); // BL

      weights.at<cv::Vec2f>(v_min, u_min) += cv::Vec2f{w0,w0};
      weights.at<cv::Vec2f>(v_min, u_max) += cv::Vec2f{w1,w1};
      weights.at<cv::Vec2f>(v_max, u_min) += cv::Vec2f{w3,w3};
      weights.at<cv::Vec2f>(v_max, u_max) += cv::Vec2f{w2,w2};

      inverse_flow.at<cv::Vec2f>(v_min, u_min) += w0 * cv::Vec2f{-du,-dv};
      inverse_flow.at<cv::Vec2f>(v_min, u_max) += w1 * cv::Vec2f{-du,-dv};
      inverse_flow.at<cv::Vec2f>(v_max, u_min) += w3 * cv::Vec2f{-du,-dv}; 
      inverse_flow.at<cv::Vec2f>(v_max, u_max) += w2 * cv::Vec2f{-du,-dv}; 
    }
  }

  cv::divide(inverse_flow, weights, inverse_flow);

  return inverse_flow;                      
}