Temporal propagation in DISOpticalFlow

Added an option to pass an initial approximation of optical flow in
DISOpticalFlow. Added a python sample that demonstrates the use of this
feature for temporal propagation of flow vectors.

modules/optflow/include/opencv2/optflow.hpp

@@ -167,7 +167,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;
 
@@ -255,6 +255,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
 {
@@ -323,7 +328,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;

modules/optflow/src/dis_flow.cpp

@@ -110,6 +110,9 @@ class DISOpticalFlowImpl : public DISOpticalFlow
     vector<Mat_<float> > Ux; //!< x component of the flow vectors
     vector<Mat_<float> > Uy; //!< y component of the flow vectors
 
+    vector<Mat_<float> > initial_Ux; //!< x component of the initial flow field, if one was passed as an input
+    vector<Mat_<float> > initial_Uy; //!< y component of the initial flow field, if one was passed as an input
+
     Mat_<Vec2f> U; //!< a buffer for the merged flow
 
     Mat_<float> Sx; //!< intermediate sparse flow representation (x component)
@@ -121,8 +124,8 @@ class DISOpticalFlowImpl : public DISOpticalFlow
     Mat_<float> I0xy_buf; //!< sum of x and y gradient products
 
     /* Extra buffers that are useful if patch mean-normalization is used: */
-    Mat_<float> I0x_buf; //!< sum of of x gradient values
-    Mat_<float> I0y_buf; //!< sum of of y gradient values
+    Mat_<float> I0x_buf; //!< sum of x gradient values
+    Mat_<float> I0y_buf; //!< sum of y gradient values
 
     /* Auxiliary buffers used in structure tensor computation: */
     Mat_<float> I0xx_buf_aux;
@@ -134,7 +137,7 @@ class DISOpticalFlowImpl : public DISOpticalFlow
     vector<Ptr<VariationalRefinement> > 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;n<range.end;n++)
-            (*this)(Range(n,n+1));
+        for (int n = range.start; n < range.end; n++)
+            (*this)(Range(n, n + 1));
         return;
     }
     int psz = dis->patch_size;
@@ -708,6 +729,15 @@ void DISOpticalFlowImpl::PatchInverseSearch_ParBody::operator()(const Range &ran
     float *x_ptr = dis->I0x_buf.ptr<float>();
     float *y_ptr = dis->I0y_buf.ptr<float>();
 
+    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<float>();
+        initial_Uy_ptr = dis->initial_Uy[pyr_level].ptr<float>();
+        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),

samples/python2/dis_opt_flow.py

+#!/usr/bin/env python
+
+'''
+example to show optical flow estimation using DISOpticalFlow
+
+USAGE: dis_opt_flow.py [<video_source>]
+
+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()