diff --git a/modules/optflow/include/opencv2/optflow.hpp b/modules/optflow/include/opencv2/optflow.hpp index e95124336..ccd7f8865 100644 --- a/modules/optflow/include/opencv2/optflow.hpp +++ b/modules/optflow/include/opencv2/optflow.hpp @@ -170,7 +170,7 @@ procedure can be found in @cite Brox2004 class CV_EXPORTS_W VariationalRefinement : public DenseOpticalFlow { public: - /** @brief calc function overload to handle separate horizontal (u) and vertical (v) flow components + /** @brief @ref calc function overload to handle separate horizontal (u) and vertical (v) flow components (to avoid extra splits/merges) */ CV_WRAP virtual void calcUV(InputArray I0, InputArray I1, InputOutputArray flow_u, InputOutputArray flow_v) = 0; @@ -258,6 +258,11 @@ This class implements the Dense Inverse Search (DIS) optical flow algorithm. Mor details about the algorithm can be found at @cite Kroeger2016 . Includes three presets with preselected parameters to provide reasonable trade-off between speed and quality. However, even the slowest preset is still relatively fast, use DeepFlow if you need better quality and don't care about speed. + +This implementation includes several additional features compared to the algorithm described in the paper, +including spatial propagation of flow vectors (@ref getUseSpatialPropagation), as well as an option to +utilize an initial flow approximation passed to @ref calc (which is, essentially, temporal propagation, +if the previous frame's flow field is passed). */ class CV_EXPORTS_W DISOpticalFlow : public DenseOpticalFlow { @@ -326,7 +331,7 @@ public: /** @brief Whether to use mean-normalization of patches when computing patch distance. It is turned on by default as it typically provides a noticeable quality boost because of increased robustness to - illumanition variations. Turn it off if you are certain that your sequence does't contain any changes + illumination variations. Turn it off if you are certain that your sequence doesn't contain any changes in illumination. @see setUseMeanNormalization */ CV_WRAP virtual bool getUseMeanNormalization() const = 0; diff --git a/modules/optflow/src/dis_flow.cpp b/modules/optflow/src/dis_flow.cpp index 15b3133e9..a485013cb 100644 --- a/modules/optflow/src/dis_flow.cpp +++ b/modules/optflow/src/dis_flow.cpp @@ -110,6 +110,9 @@ class DISOpticalFlowImpl : public DISOpticalFlow vector > Ux; //!< x component of the flow vectors vector > Uy; //!< y component of the flow vectors + vector > initial_Ux; //!< x component of the initial flow field, if one was passed as an input + vector > initial_Uy; //!< y component of the initial flow field, if one was passed as an input + Mat_ U; //!< a buffer for the merged flow Mat_ Sx; //!< intermediate sparse flow representation (x component) @@ -121,8 +124,8 @@ class DISOpticalFlowImpl : public DISOpticalFlow Mat_ I0xy_buf; //!< sum of x and y gradient products /* Extra buffers that are useful if patch mean-normalization is used: */ - Mat_ I0x_buf; //!< sum of of x gradient values - Mat_ I0y_buf; //!< sum of of y gradient values + Mat_ I0x_buf; //!< sum of x gradient values + Mat_ I0y_buf; //!< sum of y gradient values /* Auxiliary buffers used in structure tensor computation: */ Mat_ I0xx_buf_aux; @@ -134,7 +137,7 @@ class DISOpticalFlowImpl : public DISOpticalFlow vector > variational_refinement_processors; private: //!< private methods and parallel sections - void prepareBuffers(Mat &I0, Mat &I1); + void prepareBuffers(Mat &I0, Mat &I1, Mat &flow, bool use_flow); void precomputeStructureTensor(Mat &dst_I0xx, Mat &dst_I0yy, Mat &dst_I0xy, Mat &dst_I0x, Mat &dst_I0y, Mat &I0x, Mat &I0y); @@ -144,10 +147,11 @@ class DISOpticalFlowImpl : public DISOpticalFlow int nstripes, stripe_sz; int hs; Mat *Sx, *Sy, *Ux, *Uy, *I0, *I1, *I0x, *I0y; - int num_iter; + int num_iter, pyr_level; PatchInverseSearch_ParBody(DISOpticalFlowImpl &_dis, int _nstripes, int _hs, Mat &dst_Sx, Mat &dst_Sy, - Mat &src_Ux, Mat &src_Uy, Mat &_I0, Mat &_I1, Mat &_I0x, Mat &_I0y, int _num_iter); + Mat &src_Ux, Mat &src_Uy, Mat &_I0, Mat &_I1, Mat &_I0x, Mat &_I0y, int _num_iter, + int _pyr_level); void operator()(const Range &range) const; }; @@ -185,7 +189,7 @@ DISOpticalFlowImpl::DISOpticalFlowImpl() variational_refinement_processors.push_back(createVariationalFlowRefinement()); } -void DISOpticalFlowImpl::prepareBuffers(Mat &I0, Mat &I1) +void DISOpticalFlowImpl::prepareBuffers(Mat &I0, Mat &I1, Mat &flow, bool use_flow) { I0s.resize(coarsest_scale + 1); I1s.resize(coarsest_scale + 1); @@ -195,6 +199,14 @@ void DISOpticalFlowImpl::prepareBuffers(Mat &I0, Mat &I1) Ux.resize(coarsest_scale + 1); Uy.resize(coarsest_scale + 1); + Mat flow_uv[2]; + if (use_flow) + { + split(flow, flow_uv); + initial_Ux.resize(coarsest_scale + 1); + initial_Uy.resize(coarsest_scale + 1); + } + int fraction = 1; int cur_rows = 0, cur_cols = 0; @@ -237,8 +249,6 @@ void DISOpticalFlowImpl::prepareBuffers(Mat &I0, Mat &I1) resize(I1s[i - 1], I1s[i], I1s[i].size(), 0.0, 0.0, INTER_AREA); } - fraction *= 2; - if (i >= finest_scale) { I1s_ext[i].create(cur_rows + 2 * border_size, cur_cols + 2 * border_size); @@ -253,7 +263,17 @@ void DISOpticalFlowImpl::prepareBuffers(Mat &I0, Mat &I1) variational_refinement_processors[i]->setGamma(variational_refinement_gamma); variational_refinement_processors[i]->setSorIterations(5); variational_refinement_processors[i]->setFixedPointIterations(variational_refinement_iter); + + if (use_flow) + { + resize(flow_uv[0], initial_Ux[i], Size(cur_cols, cur_rows)); + initial_Ux[i] /= fraction; + resize(flow_uv[1], initial_Uy[i], Size(cur_cols, cur_rows)); + initial_Uy[i] /= fraction; + } } + + fraction *= 2; } } @@ -377,9 +397,10 @@ void DISOpticalFlowImpl::precomputeStructureTensor(Mat &dst_I0xx, Mat &dst_I0yy, DISOpticalFlowImpl::PatchInverseSearch_ParBody::PatchInverseSearch_ParBody(DISOpticalFlowImpl &_dis, int _nstripes, int _hs, Mat &dst_Sx, Mat &dst_Sy, Mat &src_Ux, Mat &src_Uy, Mat &_I0, Mat &_I1, - Mat &_I0x, Mat &_I0y, int _num_iter) + Mat &_I0x, Mat &_I0y, int _num_iter, + int _pyr_level) : dis(&_dis), nstripes(_nstripes), hs(_hs), Sx(&dst_Sx), Sy(&dst_Sy), Ux(&src_Ux), Uy(&src_Uy), I0(&_I0), I1(&_I1), - I0x(&_I0x), I0y(&_I0y), num_iter(_num_iter) + I0x(&_I0x), I0y(&_I0y), num_iter(_num_iter), pyr_level(_pyr_level) { stripe_sz = (int)ceil(hs / (double)nstripes); } @@ -676,10 +697,10 @@ inline float computeSSDMeanNorm(uchar *I0_ptr, uchar *I1_ptr, int I0_stride, int void DISOpticalFlowImpl::PatchInverseSearch_ParBody::operator()(const Range &range) const { // force separate processing of stripes if we are using spatial propagation: - if(dis->use_spatial_propagation && range.end>range.start+1) + if (dis->use_spatial_propagation && range.end > range.start + 1) { - for(int n=range.start;npatch_size; @@ -708,6 +729,15 @@ void DISOpticalFlowImpl::PatchInverseSearch_ParBody::operator()(const Range &ran float *x_ptr = dis->I0x_buf.ptr(); float *y_ptr = dis->I0y_buf.ptr(); + bool use_temporal_candidates = false; + float *initial_Ux_ptr = NULL, *initial_Uy_ptr = NULL; + if (!dis->initial_Ux.empty()) + { + initial_Ux_ptr = dis->initial_Ux[pyr_level].ptr(); + initial_Uy_ptr = dis->initial_Uy[pyr_level].ptr(); + use_temporal_candidates = true; + } + int i, j, dir; int start_is, end_is, start_js, end_js; int start_i, start_j; @@ -772,11 +802,28 @@ void DISOpticalFlowImpl::PatchInverseSearch_ParBody::operator()(const Range &ran Sy_ptr[is * dis->ws + js] = Uy_ptr[(i + psz2) * dis->w + j + psz2]; } - if (dis->use_spatial_propagation) + float min_SSD = INF, cur_SSD; + if (use_temporal_candidates || dis->use_spatial_propagation) { - /* Updating the current Sx_ptr, Sy_ptr to the best candidate: */ - float min_SSD, cur_SSD; COMPUTE_SSD(min_SSD, Sx_ptr[is * dis->ws + js], Sy_ptr[is * dis->ws + js]); + } + + if (use_temporal_candidates) + { + /* Try temporal candidates (vectors from the initial flow field that was passed to the function) */ + COMPUTE_SSD(cur_SSD, initial_Ux_ptr[(i + psz2) * dis->w + j + psz2], + initial_Uy_ptr[(i + psz2) * dis->w + j + psz2]); + if (cur_SSD < min_SSD) + { + min_SSD = cur_SSD; + Sx_ptr[is * dis->ws + js] = initial_Ux_ptr[(i + psz2) * dis->w + j + psz2]; + Sy_ptr[is * dis->ws + js] = initial_Uy_ptr[(i + psz2) * dis->w + j + psz2]; + } + } + + if (dis->use_spatial_propagation) + { + /* Try spatial candidates: */ if (dir * js > dir * start_js) { COMPUTE_SSD(cur_SSD, Sx_ptr[is * dis->ws + js - dir], Sy_ptr[is * dis->ws + js - dir]); @@ -967,12 +1014,16 @@ void DISOpticalFlowImpl::calc(InputArray I0, InputArray I1, InputOutputArray flo Mat I0Mat = I0.getMat(); Mat I1Mat = I1.getMat(); - flow.create(I1Mat.size(), CV_32FC2); + bool use_input_flow = false; + if (flow.sameSize(I0) && flow.depth() == CV_32F && flow.channels() == 2) + use_input_flow = true; + else + flow.create(I1Mat.size(), CV_32FC2); Mat &flowMat = flow.getMatRef(); coarsest_scale = (int)(log((2 * I0Mat.cols) / (4.0 * patch_size)) / log(2.0) + 0.5) - 1; int num_stripes = getNumThreads(); - prepareBuffers(I0Mat, I1Mat); + prepareBuffers(I0Mat, I1Mat, flowMat, use_input_flow); Ux[coarsest_scale].setTo(0.0f); Uy[coarsest_scale].setTo(0.0f); @@ -990,13 +1041,13 @@ void DISOpticalFlowImpl::calc(InputArray I0, InputArray I1, InputOutputArray flo * with spatial propagation reproducible */ parallel_for_(Range(0, 8), PatchInverseSearch_ParBody(*this, 8, hs, Sx, Sy, Ux[i], Uy[i], I0s[i], - I1s_ext[i], I0xs[i], I0ys[i], 2)); + I1s_ext[i], I0xs[i], I0ys[i], 2, i)); } else { parallel_for_(Range(0, num_stripes), PatchInverseSearch_ParBody(*this, num_stripes, hs, Sx, Sy, Ux[i], Uy[i], I0s[i], I1s_ext[i], - I0xs[i], I0ys[i], 1)); + I0xs[i], I0ys[i], 1, i)); } parallel_for_(Range(0, num_stripes), diff --git a/samples/python2/dis_opt_flow.py b/samples/python2/dis_opt_flow.py new file mode 100644 index 000000000..731a3aa5b --- /dev/null +++ b/samples/python2/dis_opt_flow.py @@ -0,0 +1,114 @@ +#!/usr/bin/env python + +''' +example to show optical flow estimation using DISOpticalFlow + +USAGE: dis_opt_flow.py [] + +Keys: + 1 - toggle HSV flow visualization + 2 - toggle glitch + 3 - toggle spatial propagation of flow vectors + 4 - toggle temporal propagation of flow vectors +ESC - exit +''' + +# Python 2/3 compatibility +from __future__ import print_function + +import numpy as np +import cv2 +import video + + +def draw_flow(img, flow, step=16): + h, w = img.shape[:2] + y, x = np.mgrid[step/2:h:step, step/2:w:step].reshape(2,-1).astype(int) + fx, fy = flow[y,x].T + lines = np.vstack([x, y, x+fx, y+fy]).T.reshape(-1, 2, 2) + lines = np.int32(lines + 0.5) + vis = cv2.cvtColor(img, cv2.COLOR_GRAY2BGR) + cv2.polylines(vis, lines, 0, (0, 255, 0)) + for (x1, y1), (x2, y2) in lines: + cv2.circle(vis, (x1, y1), 1, (0, 255, 0), -1) + return vis + + +def draw_hsv(flow): + h, w = flow.shape[:2] + fx, fy = flow[:,:,0], flow[:,:,1] + ang = np.arctan2(fy, fx) + np.pi + v = np.sqrt(fx*fx+fy*fy) + hsv = np.zeros((h, w, 3), np.uint8) + hsv[...,0] = ang*(180/np.pi/2) + hsv[...,1] = 255 + hsv[...,2] = np.minimum(v*4, 255) + bgr = cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR) + return bgr + + +def warp_flow(img, flow): + h, w = flow.shape[:2] + flow = -flow + flow[:,:,0] += np.arange(w) + flow[:,:,1] += np.arange(h)[:,np.newaxis] + res = cv2.remap(img, flow, None, cv2.INTER_LINEAR) + return res + + +if __name__ == '__main__': + import sys + print(__doc__) + try: + fn = sys.argv[1] + except IndexError: + fn = 0 + + cam = video.create_capture(fn) + ret, prev = cam.read() + prevgray = cv2.cvtColor(prev, cv2.COLOR_BGR2GRAY) + show_hsv = False + show_glitch = False + use_spatial_propagation = False + use_temporal_propagation = True + cur_glitch = prev.copy() + inst = cv2.optflow.createOptFlow_DIS(cv2.optflow.DISOPTICAL_FLOW_PRESET_MEDIUM) + inst.setUseSpatialPropagation(use_spatial_propagation) + + flow = None + while True: + ret, img = cam.read() + gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) + if flow is not None and use_temporal_propagation: + #warp previous flow to get an initial approximation for the current flow: + flow = inst.calc(prevgray, gray, warp_flow(flow,flow)) + else: + flow = inst.calc(prevgray, gray, None) + prevgray = gray + + cv2.imshow('flow', draw_flow(gray, flow)) + if show_hsv: + cv2.imshow('flow HSV', draw_hsv(flow)) + if show_glitch: + cur_glitch = warp_flow(cur_glitch, flow) + cv2.imshow('glitch', cur_glitch) + + ch = 0xFF & cv2.waitKey(5) + if ch == 27: + break + if ch == ord('1'): + show_hsv = not show_hsv + print('HSV flow visualization is', ['off', 'on'][show_hsv]) + if ch == ord('2'): + show_glitch = not show_glitch + if show_glitch: + cur_glitch = img.copy() + print('glitch is', ['off', 'on'][show_glitch]) + if ch == ord('3'): + use_spatial_propagation = not use_spatial_propagation + inst.setUseSpatialPropagation(use_spatial_propagation) + print('spatial propagation is', ['off', 'on'][use_spatial_propagation]) + if ch == ord('4'): + use_temporal_propagation = not use_temporal_propagation + print('temporal propagation is', ['off', 'on'][use_temporal_propagation]) + cv2.destroyAllWindows()