"SimpleFlow" optical flow estimation algorithm (GSoC project) declaration in…

"SimpleFlow" optical flow estimation algorithm (GSoC project) declaration in includes, implementation, usage example, test

modules/video/doc/motion_analysis_and_object_tracking.rst

@@ -597,6 +597,48 @@ Returns background image
 See :ocv:func:`BackgroundSubtractor::getBackgroundImage`.
 
 
+calcOpticalFlowSF
+-----------
+Calculate an optical flow using "SimpleFlow" algorithm.
+
+.. ocv:function:: void calcOpticalFlowSF( Mat& prev, Mat& next, Mat& flowX, Mat& flowY, int layers, int averaging_block_size, int max_flow, double sigma_dist, double sigma_color, int postprocess_window, double sigma_dist_fix, double sigma_color_fix, double occ_thr, int upscale_averaging_radiud, double upscale_sigma_dist, double upscale_sigma_color, double speed_up_thr)
+
+    :param prev: First 8-bit 3-channel image.
+
+    :param next: Second 8-bit 3-channel image 
+
+    :param flowX: X-coordinate of estimated flow
+
+    :param flowY: Y-coordinate of estimated flow
+
+    :param layers: Number of layers
+
+    :param averaging_block_size: Size of block through which we sum up when calculate cost function for pixel
+
+    :param max_flow: maximal flow that we search at each level
+
+    :param sigma_dist: vector smooth spatial sigma parameter
+
+    :param sigma_color: vector smooth color sigma parameter
+
+    :param postprocess_window: window size for postprocess cross bilateral filter
+
+    :param sigma_dist_fix: spatial sigma for postprocess cross bilateralf filter
+
+    :param sigma_color_fix: color sigma for postprocess cross bilateral filter
+
+    :param occ_thr: threshold for detecting occlusions
+
+    :param upscale_averaging_radiud: window size for bilateral upscale operation
+
+    :param upscale_sigma_dist: spatial sigma for bilateral upscale operation
+
+    :param upscale_sigma_color: color sigma for bilateral upscale operation
+
+    :param speed_up_thr: threshold to detect point with irregular flow - where flow should be recalculated after upscale
+
+See [Tao2012]_. And site of project - http://graphics.berkeley.edu/papers/Tao-SAN-2012-05/. 
+
 .. [Bouguet00] Jean-Yves Bouguet. Pyramidal Implementation of the Lucas Kanade Feature Tracker.
 
 .. [Bradski98] Bradski, G.R. "Computer Vision Face Tracking for Use in a Perceptual User Interface", Intel, 1998
@@ -612,3 +654,5 @@ See :ocv:func:`BackgroundSubtractor::getBackgroundImage`.
 .. [Lucas81] Lucas, B., and Kanade, T. An Iterative Image Registration Technique with an Application to Stereo Vision, Proc. of 7th International Joint Conference on Artificial Intelligence (IJCAI), pp. 674-679.
 
 .. [Welch95] Greg Welch and Gary Bishop “An Introduction to the Kalman Filter”, 1995
+
+.. [Tao2012] Michael Tao, Jiamin Bai, Pushmeet Kohli and Sylvain Paris. SimpleFlow: A Non-iterative, Sublinear Optical Flow Algorithm. Computer Graphics Forum (Eurographics 2012)

modules/video/include/opencv2/video/tracking.hpp

@@ -326,7 +326,26 @@ CV_EXPORTS_W void calcOpticalFlowFarneback( InputArray prev, InputArray next,
 // that maps one 2D point set to another or one image to another.
 CV_EXPORTS_W Mat estimateRigidTransform( InputArray src, InputArray dst,
                                          bool fullAffine);
-    
+  
+//! computes dense optical flow using Simple Flow algorithm
+CV_EXPORTS_W void calcOpticalFlowSF(Mat& from, 
+                                    Mat& to, 
+                                    Mat& flowX,
+                                    Mat& flowY,
+                                    int layers,
+                                    int averaging_block_size, 
+                                    int max_flow,
+                                    double sigma_dist,
+                                    double sigma_color,
+                                    int postprocess_window, 
+                                    double sigma_dist_fix, 
+                                    double sigma_color_fix,
+                                    double occ_thr,
+                                    int upscale_averaging_radius,
+                                    double upscale_sigma_dist,
+                                    double upscale_sigma_color,
+                                    double speed_up_thr);
+
 }
 
 #endif

modules/video/src/simpleflow.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                          License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_SIMPLEFLOW_H__
+#define __OPENCV_SIMPLEFLOW_H__
+
+#include <vector>
+
+using namespace std;
+
+#define MASK_TRUE_VALUE 255
+#define UNKNOWN_FLOW_THRESH 1e9
+
+namespace cv {
+
+struct Flow {
+  Mat u, v;
+
+  Flow() {;}
+
+  Flow(Mat& _u, Mat& _v)
+    : u(_u), v(_v) {;}
+  
+  Flow(int rows, int cols) {
+    u = Mat::zeros(rows, cols, CV_64F);
+    v = Mat::zeros(rows, cols, CV_64F);
+  }
+};
+
+inline static double dist(const Vec3b& p1, const Vec3b& p2) {
+  return (p1[0] - p2[0]) * (p1[0] - p2[0]) +
+         (p1[1] - p2[1]) * (p1[1] - p2[1]) +
+         (p1[2] - p2[2]) * (p1[2] - p2[2]);
+}
+
+inline static double dist(const Point2f& p1, const Point2f& p2) {
+  return (p1.x - p2.x) * (p1.x - p2.x) +
+         (p1.y - p2.y) * (p1.y - p2.y);
+}
+
+inline static double dist(double x1, double y1, double x2, double y2) {
+  return (x1 - x2) * (x1 - x2) +
+         (y1 - y2) * (y1 - y2);
+}
+
+inline static int dist(int x1, int y1, int x2, int y2) {
+  return (x1 - x2) * (x1 - x2) +
+         (y1 - y2) * (y1 - y2);
+}
+
+template<class T>
+inline static T min(T t1, T t2, T t3) {
+  return (t1 <= t2 && t1 <= t3) ? t1 : min(t2, t3);
+}
+
+template<class T>
+vector<vector<T> > build(int n, int m) {
+  vector<vector<T> > res(n);
+  for (int i = 0; i < n; ++i) {
+    res[i].resize(m, 0);
+  }
+  return res;
+}
+
+class WeightedCrossBilateralFilter {
+public:
+  WeightedCrossBilateralFilter(const Mat& _image,
+                       int _windowSize,
+                       double _sigmaDist,
+                       double _sigmaColor);
+
+  Mat apply(Mat& matrix, Mat& weights);
+
+private:
+  double convolution(Mat& matrix, int row, int col, Mat& weights);
+
+  Mat image;
+  int windowSize;
+  double sigmaDist, sigmaColor;
+
+  vector<double> expDist;
+  vector<vector<vector<vector<double> > > > wc;
+};
+}
+
+#endif

modules/video/test/test_simpleflow.cpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "test_precomp.hpp"
+
+#include <string>
+
+using namespace std;
+
+/* ///////////////////// simpleflow_test ///////////////////////// */
+
+class CV_SimpleFlowTest : public cvtest::BaseTest
+{
+public:
+    CV_SimpleFlowTest();
+protected:
+    void run(int);
+};
+
+
+CV_SimpleFlowTest::CV_SimpleFlowTest() {}
+
+static void readOpticalFlowFromFile(FILE* file, cv::Mat& flowX, cv::Mat& flowY) {
+  char header[5];
+  if (fread(header, 1, 4, file) < 4 && (string)header != "PIEH") {
+    return;
+  }
+
+  int cols, rows;
+  if (fread(&cols, sizeof(int), 1, file) != 1||
+      fread(&rows, sizeof(int), 1, file) != 1) {
+    return;
+  }
+
+  flowX = cv::Mat::zeros(rows, cols, CV_64F);
+  flowY = cv::Mat::zeros(rows, cols, CV_64F);
+
+  for (int i = 0; i < rows; ++i) {
+    for (int j = 0; j < cols; ++j) {
+      float uPoint, vPoint;
+      if (fread(&uPoint, sizeof(float), 1, file) != 1 ||
+          fread(&vPoint, sizeof(float), 1, file) != 1) {
+        flowX.release();
+        flowY.release();
+        return;
+      }
+  
+      flowX.at<double>(i, j) = uPoint;
+      flowY.at<double>(i, j) = vPoint;
+    }
+  }
+}
+
+static bool isFlowCorrect(double u) {
+  return !isnan(u) && (fabs(u) < 1e9);
+}
+
+static double calc_rmse(cv::Mat flow1X, cv::Mat flow1Y, cv::Mat flow2X, cv::Mat flow2Y) {
+  long double sum;
+  int counter = 0;
+  const int rows = flow1X.rows;
+  const int cols = flow1X.cols;
+
+  for (int y = 0; y < rows; ++y) {
+    for (int x = 0; x < cols; ++x) {
+      double u1 = flow1X.at<double>(y, x);
+      double v1 = flow1Y.at<double>(y, x);
+      double u2 = flow2X.at<double>(y, x);
+      double v2 = flow2Y.at<double>(y, x);
+      if (isFlowCorrect(u1) && isFlowCorrect(u2) && isFlowCorrect(v1) && isFlowCorrect(v2)) {
+        sum += (u1-u2)*(u1-u2) + (v1-v2)*(v1-v2);
+        counter++;
+      }
+    }
+  }
+  return sqrt((double)sum / (1e-9 + counter));
+}
+
+void CV_SimpleFlowTest::run(int) {
+    int code = cvtest::TS::OK;
+    
+    const double MAX_RMSE = 0.6;
+    const string frame1_path = ts->get_data_path() + "optflow/RubberWhale1.png";
+    const string frame2_path = ts->get_data_path() + "optflow/RubberWhale2.png";
+    const string gt_flow_path = ts->get_data_path() + "optflow/RubberWhale.flo";
+
+    cv::Mat frame1 = cv::imread(frame1_path);
+    cv::Mat frame2 = cv::imread(frame2_path);
+
+    if (frame1.empty()) {
+      ts->printf(cvtest::TS::LOG, "could not read image %s\n", frame2_path.c_str());
+      ts->set_failed_test_info(cvtest::TS::FAIL_MISSING_TEST_DATA);
+      return;
+    }
+    
+    if (frame2.empty()) {
+      ts->printf(cvtest::TS::LOG, "could not read image %s\n", frame2_path.c_str());
+      ts->set_failed_test_info(cvtest::TS::FAIL_MISSING_TEST_DATA);
+      return;
+    }
+
+    if (frame1.rows != frame2.rows && frame1.cols != frame2.cols) {
+      ts->printf(cvtest::TS::LOG, "images should be of equal sizes (%s and %s)",
+                 frame1_path.c_str(), frame2_path.c_str());
+      ts->set_failed_test_info(cvtest::TS::FAIL_MISSING_TEST_DATA);
+      return;
+    }
+
+    if (frame1.type() != 16 || frame2.type() != 16) {
+      ts->printf(cvtest::TS::LOG, "images should be of equal type CV_8UC3 (%s and %s)",
+                 frame1_path.c_str(), frame2_path.c_str());
+      ts->set_failed_test_info(cvtest::TS::FAIL_MISSING_TEST_DATA);
+      return;
+    }
+
+    cv::Mat flowX_gt, flowY_gt;
+
+    FILE* gt_flow_file = fopen(gt_flow_path.c_str(), "rb");
+    if (gt_flow_file == NULL) {
+      ts->printf(cvtest::TS::LOG, "could not read ground-thuth flow from file %s",
+                 gt_flow_path.c_str());
+      ts->set_failed_test_info(cvtest::TS::FAIL_MISSING_TEST_DATA);
+      return;
+    }
+    readOpticalFlowFromFile(gt_flow_file, flowX_gt, flowY_gt);
+    if (flowX_gt.empty() || flowY_gt.empty()) {
+      ts->printf(cvtest::TS::LOG, "error while reading flow data from file %s",
+                 gt_flow_path.c_str());
+      ts->set_failed_test_info(cvtest::TS::FAIL_MISSING_TEST_DATA);
+      return;
+    }
+    fclose(gt_flow_file);
+
+    cv::Mat flowX, flowY;
+    cv::calcOpticalFlowSF(frame1, frame2, 
+                          flowX, flowY,
+                          3, 4, 2, 4.1, 25.5, 18, 55.0, 25.5, 0.35, 18, 55.0, 25.5, 10);
+
+    double rmse = calc_rmse(flowX_gt, flowY_gt, flowX, flowY);
+    
+    ts->printf(cvtest::TS::LOG, "Optical flow estimation RMSE for SimpleFlow algorithm : %lf\n",
+               rmse);
+
+    if (rmse > MAX_RMSE) {
+      ts->printf( cvtest::TS::LOG,
+                 "Too big rmse error : %lf ( >= %lf )\n", rmse, MAX_RMSE);
+      ts->set_failed_test_info(cvtest::TS::FAIL_BAD_ACCURACY);
+      return;
+    }
+}
+
+
+TEST(Video_OpticalFlowSimpleFlow, accuracy) { CV_SimpleFlowTest test; test.safe_run(); }
+
+/* End of file. */

samples/cpp/simpleflow_demo.cpp

+#include "opencv2/video/tracking.hpp"
+#include "opencv2/imgproc/imgproc.hpp"
+#include "opencv2/highgui/highgui.hpp"
+
+#include <cstdio>
+#include <iostream>
+
+using namespace cv;
+using namespace std;
+
+static void help()
+{
+    // print a welcome message, and the OpenCV version
+    printf("This is a demo of SimpleFlow optical flow algorithm,\n"
+           "Using OpenCV version %s\n\n", CV_VERSION);
+
+    printf("Usage: simpleflow_demo frame1 frame2 output_flow"
+           "\nApplication will write estimated flow "
+           "\nbetween 'frame1' and 'frame2' in binary format"
+           "\ninto file 'output_flow'"
+           "\nThen one can use code from http://vision.middlebury.edu/flow/data/"
+           "\nto convert flow in binary file to image\n");
+}
+
+// binary file format for flow data specified here:
+// http://vision.middlebury.edu/flow/data/
+static void writeOpticalFlowToFile(const Mat& u, const Mat& v, FILE* file) {
+  int cols = u.cols;
+  int rows = u.rows;
+  
+  fprintf(file, "PIEH");
+   
+  if (fwrite(&cols, sizeof(int), 1, file) != 1 ||
+      fwrite(&rows, sizeof(int), 1, file) != 1) {
+    fprintf(stderr, "writeOpticalFlowToFile : problem writing header\n");
+    exit(1);
+  }
+
+  for (int i= 0; i < u.rows; ++i) {
+    for (int j = 0; j < u.cols; ++j) {
+      float uPoint = u.at<double>(i, j);
+      float vPoint = v.at<double>(i, j);
+
+      if (fwrite(&uPoint, sizeof(float), 1, file) != 1 ||
+          fwrite(&vPoint, sizeof(float), 1, file) != 1) {
+        fprintf(stderr, "writeOpticalFlowToFile : problem writing data\n");
+        exit(1);
+      }
+    }
+  }
+}
+int main(int argc, char** argv) {
+    help();
+
+    if (argc < 4) {
+      fprintf(stderr, "Wrong number of command line arguments : %d (expected %d)\n", argc, 4);
+      exit(1);
+    }
+    
+    Mat frame1 = imread(argv[1]);
+    Mat frame2 = imread(argv[2]);
+
+    if (frame1.empty() || frame2.empty()) {
+      fprintf(stderr, "simpleflow_demo : Images cannot be read\n");
+      exit(1);
+    }
+
+    if (frame1.rows != frame2.rows && frame1.cols != frame2.cols) {
+      fprintf(stderr, "simpleflow_demo : Images should be of equal sizes\n");
+      exit(1);
+    }
+
+    if (frame1.type() != 16 || frame2.type() != 16) {
+      fprintf(stderr, "simpleflow_demo : Images should be of equal type CV_8UC3\n");
+      exit(1);
+    }
+
+    printf("simpleflow_demo : Read two images of size [rows = %d, cols = %d]\n", 
+           frame1.rows, frame1.cols);
+
+    Mat flowX, flowY;
+
+    calcOpticalFlowSF(frame1, frame2, 
+                      flowX, flowY,
+                      3, 2, 4, 4.1, 25.5, 18, 55.0, 25.5, 0.35, 18, 55.0, 25.5, 10);
+
+  FILE* file = fopen(argv[3], "wb");
+  if (file == NULL) {
+    fprintf(stderr, "simpleflow_demo : Unable to open file '%s' for writing\n", argv[3]); 
+    exit(1);
+  }
+  printf("simpleflow_demo : Writing to file\n");
+  writeOpticalFlowToFile(flowX, flowY, file);
+  fclose(file);
+  return 0;
+}