Merge pull request #68 from tpietruszka/deepflow

Deepflow implementation (and optical flow testing again)

modules/optflow/doc/dense_optflow.rst

@@ -48,6 +48,69 @@ See [Tao2012]_. And site of project - http://graphics.berkeley.edu/papers/Tao-SA
 .. note::
 
    * An example using the simpleFlow algorithm can be found at samples/simpleflow_demo.cpp
-   
+
+
+optflow::OpticalFlowDeepFlow
+-----------------------------
+DeepFlow optical flow algorithm implementation.
+
+.. ocv:class:: optflow::OpticalFlowDeepFlow: public DenseOpticalFlow
+
+
+  The class implements the DeepFlow optical flow algorithm described in [Weinzaepfel2013]_ . See also http://lear.inrialpes.fr/src/deepmatching/ .
+
+  Parameters - class fields - that may be modified after creating a class instance:
+
+    .. ocv:member:: float alpha
+
+    Smoothness assumption weight
+
+    .. ocv:member:: float delta
+
+    Color constancy assumption weight
+
+    .. ocv:member:: float gamma
+
+    Gradient constancy weight
+        
+    .. ocv:member:: float sigma
+
+    Gaussian smoothing parameter
+        
+    .. ocv:member:: int minSize
+
+    Minimal dimension of an image in the pyramid (next, smaller images in the pyramid are
+    generated until one of the dimensions reaches this size)
+        
+    .. ocv:member:: float downscaleFactor
+
+    Scaling factor in the image pyramid (must be < 1) 
+        
+    .. ocv:member:: int fixedPointIterations
+
+    How many iterations on each level of the pyramid
+        
+    .. ocv:member:: int sorIterations
+
+    Iterations of Succesive Over-Relaxation (solver)
+        
+    .. ocv:member:: float omega
+
+    Relaxation factor in SOR
+
+
+optflow::createOptFlow_DeepFlow
+---------------------------------
+Create an OpticalFlowDeepFlow instance
+
+  .. ocv:function:: Ptr<DenseOpticalFlow> optflow::createOptFlow_DeepFlow()
+
+  Returns a pointer to a ``DenseOpticalFlow`` instance - see ``DenseOpticalFlow::calc()``
+
+
+  
+
+        
 
 .. [Tao2012] Michael Tao, Jiamin Bai, Pushmeet Kohli and Sylvain Paris. SimpleFlow: A Non-iterative, Sublinear Optical Flow Algorithm. Computer Graphics Forum (Eurographics 2012)
+.. [Weinzaepfel2013] P. Weinzaepfel, J. Revaud, Z. Harchaoui, and C. Schmid, “DeepFlow: Large Displacement Optical Flow with Deep Matching,” 2013 IEEE Int. Conf. Comput. Vis., pp. 1385–1392, Dec. 2013.

modules/optflow/doc/optflow.rst

@@ -11,3 +11,4 @@ The opencv_optflow module includes different algorithms for tracking points
 
     dense_optflow
     motion_templates
+    optflow_io

modules/optflow/doc/optflow_io.rst

+Optical flow input / output
+============================
+
+Functions reading and writing .flo files in "Middlebury" format, see: [MiddleburyFlo]_ 
+
+readOpticalFlow
+-----------------
+Read a .flo file
+
+.. ocv:function:: Mat readOpticalFlow( const String& path )
+    
+    :param path: Path to the file to be loaded
+    
+The function readOpticalFlow loads a flow field from a file and returns it as a single matrix.
+Resulting ``Mat`` has a type ``CV_32FC2`` - floating-point, 2-channel. 
+First channel corresponds to the flow in the horizontal direction (u), second - vertical (v).
+
+writeOpticalFlow
+-----------------
+Write a .flo to disk
+
+.. ocv:function:: bool writeOpticalFlow( const String& path, InputArray flow )
+
+    :param path: Path to the file to be written
+    
+    :param flow: Flow field to be stored
+    
+ The function stores a flow field in a file, returns ``true`` on success, ``false`` otherwise.
+
+ The flow field must be a 2-channel, floating-point matrix (``CV_32FC2``). 
+ First channel corresponds to the flow in the horizontal direction (u), second - vertical (v).
+ 
+ 
+.. [MiddleburyFlo] http://vision.middlebury.edu/flow/code/flow-code/README.txt
+ 
+ 
\ No newline at end of file

modules/optflow/include/opencv2/optflow.hpp

@@ -41,6 +41,7 @@ the use of this software, even if advised of the possibility of such damage.
 #define __OPENCV_OPTFLOW_HPP__
 
 #include "opencv2/core.hpp"
+#include "opencv2/video.hpp"
 
 namespace cv
 {
@@ -57,8 +58,26 @@ CV_EXPORTS_W void calcOpticalFlowSF( InputArray from, InputArray to, OutputArray
                                      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 );    
-    
-}
+
+//! reads optical flow from a file, Middlebury format:
+// http://vision.middlebury.edu/flow/code/flow-code/README.txt
+CV_EXPORTS_W Mat readOpticalFlow( const String& path );
+//! writes optical flow to a file, Middlebury format
+CV_EXPORTS_W bool writeOpticalFlow( const String& path, InputArray flow );
+
+
+
+
+// DeepFlow implementation, based on:
+//   P. Weinzaepfel, J. Revaud, Z. Harchaoui, and C. Schmid, “DeepFlow: Large Displacement Optical Flow with Deep Matching,”
+CV_EXPORTS_W Ptr<DenseOpticalFlow> createOptFlow_DeepFlow();
+
+// Additional interface to the SimpleFlow algorithm - calcOpticalFlowSF()
+CV_EXPORTS_W Ptr<DenseOpticalFlow> createOptFlow_SimpleFlow();
+
+// Additional interface to the Farneback's algorithm - calcOpticalFlowFarneback()
+CV_EXPORTS_W Ptr<DenseOpticalFlow> createOptFlow_Farneback();
+} //optflow
 }
 
 #include "opencv2/optflow/motempl.hpp"

modules/optflow/samples/optical_flow_evaluation.cpp

+#include "opencv2/highgui.hpp"
+#include "opencv2/video.hpp"
+#include "opencv2/optflow.hpp"
+#include <fstream>
+
+using namespace std;
+using namespace cv;
+using namespace optflow;
+
+const String keys = "{help h usage ? |      | print this message   }"
+        "{@image1        |      | image1               }"
+        "{@image2        |      | image2               }"
+        "{@algorithm     |      | [farneback, simpleflow, tvl1 or deepflow] }"
+        "{@groundtruth   |      | path to the .flo file  (optional), Middlebury format }"
+        "{m measure      |endpoint| error measure - [endpoint or angular] }"
+        "{r region       |all   | region to compute stats about [all, discontinuities, untextured] }"
+        "{d display      |      | display additional info images (pauses program execution) }";
+
+inline bool isFlowCorrect( const Point2f u )
+{
+    return !cvIsNaN(u.x) && !cvIsNaN(u.y) && (fabs(u.x) < 1e9) && (fabs(u.y) < 1e9);
+}
+inline bool isFlowCorrect( const Point3f u )
+{
+    return !cvIsNaN(u.x) && !cvIsNaN(u.y) && !cvIsNaN(u.z) && (fabs(u.x) < 1e9) && (fabs(u.y) < 1e9)
+            && (fabs(u.z) < 1e9);
+}
+static Mat endpointError( const Mat_<Point2f>& flow1, const Mat_<Point2f>& flow2 )
+{
+    Mat result(flow1.size(), CV_32FC1);
+    for ( int i = 0; i < flow1.rows; ++i )
+    {
+        for ( int j = 0; j < flow1.cols; ++j )
+        {
+            const Point2f u1 = flow1(i, j);
+            const Point2f u2 = flow2(i, j);
+
+            if ( isFlowCorrect(u1) && isFlowCorrect(u2) )
+            {
+                const Point2f diff = u1 - u2;
+                result.at<float>(i, j) = sqrt((float)diff.ddot(diff)); //distance
+            } else
+                result.at<float>(i, j) = NAN;
+        }
+    }
+    return result;
+}
+static Mat angularError( const Mat_<Point2f>& flow1, const Mat_<Point2f>& flow2 )
+{
+    Mat result(flow1.size(), CV_32FC1);
+
+    for ( int i = 0; i < flow1.rows; ++i )
+    {
+        for ( int j = 0; j < flow1.cols; ++j )
+        {
+            const Point2f u1_2d = flow1(i, j);
+            const Point2f u2_2d = flow2(i, j);
+            const Point3f u1(u1_2d.x, u1_2d.y, 1);
+            const Point3f u2(u2_2d.x, u2_2d.y, 1);
+
+            if ( isFlowCorrect(u1) && isFlowCorrect(u2) )
+                result.at<float>(i, j) = acos((float)(u1.ddot(u2) / norm(u1) * norm(u2)));
+            else
+                result.at<float>(i, j) = NAN;
+        }
+    }
+    return result;
+}
+// what fraction of pixels have errors higher than given threshold?
+static float stat_RX( Mat errors, float threshold, Mat mask )
+{
+    CV_Assert(errors.size() == mask.size());
+    CV_Assert(mask.depth() == CV_8U);
+
+    int count = 0, all = 0;
+    for ( int i = 0; i < errors.rows; ++i )
+    {
+        for ( int j = 0; j < errors.cols; ++j )
+        {
+            if ( mask.at<char>(i, j) != 0 )
+            {
+                ++all;
+                if ( errors.at<float>(i, j) > threshold )
+                    ++count;
+            }
+        }
+    }
+    return (float)count / all;
+}
+static float stat_AX( Mat hist, int cutoff_count, float max_value )
+{
+    int counter = 0;
+    int bin = 0;
+    int bin_count = hist.rows;
+    while ( bin < bin_count && counter < cutoff_count )
+    {
+        counter += (int) hist.at<float>(bin, 0);
+        ++bin;
+    }
+    return (float) bin / bin_count * max_value;
+}
+static void calculateStats( Mat errors, Mat mask = Mat(), bool display_images = false )
+{
+    float R_thresholds[] = { 0.5f, 1.f, 2.f, 5.f, 10.f };
+    float A_thresholds[] = { 0.5f, 0.75f, 0.95f };
+    if ( mask.empty() )
+        mask = Mat::ones(errors.size(), CV_8U);
+    CV_Assert(errors.size() == mask.size());
+    CV_Assert(mask.depth() == CV_8U);
+
+    //displaying the mask
+    if(display_images)
+    {
+        namedWindow( "Region mask", WINDOW_AUTOSIZE );
+        imshow( "Region mask", mask );
+    }
+
+    //mean and std computation
+    Scalar s_mean, s_std;
+    float mean, std;
+    meanStdDev(errors, s_mean, s_std, mask);
+    mean = (float)s_mean[0];
+    std = (float)s_std[0];
+    printf("Average: %.2f\nStandard deviation: %.2f\n", mean, std);
+
+    //RX stats - displayed in percent
+    float R;
+    int R_thresholds_count = sizeof(R_thresholds) / sizeof(float);
+    for ( int i = 0; i < R_thresholds_count; ++i )
+    {
+        R = stat_RX(errors, R_thresholds[i], mask);
+        printf("R%.1f: %.2f%%\n", R_thresholds[i], R * 100);
+    }
+
+    //AX stats
+    double max_value;
+    minMaxLoc(errors, NULL, &max_value, NULL, NULL, mask);
+
+    Mat hist;
+    const int n_images = 1;
+    const int channels[] = { 0 };
+    const int n_dimensions = 1;
+    const int hist_bins[] = { 1024 };
+    const float iranges[] = { 0, (float) max_value };
+    const float* ranges[] = { iranges };
+    const bool uniform = true;
+    const bool accumulate = false;
+    calcHist(&errors, n_images, channels, mask, hist, n_dimensions, hist_bins, ranges, uniform,
+            accumulate);
+    int all_pixels = countNonZero(mask);
+    int cutoff_count;
+    float A;
+    int A_thresholds_count = sizeof(A_thresholds) / sizeof(float);
+    for ( int i = 0; i < A_thresholds_count; ++i )
+    {
+        cutoff_count = (int) (floor(A_thresholds[i] * all_pixels + 0.5f));
+        A = stat_AX(hist, cutoff_count, (float) max_value);
+        printf("A%.2f: %.2f\n", A_thresholds[i], A);
+    }
+}
+
+static Mat flowToDisplay(const Mat flow)
+{
+    Mat flow_split[2];
+    Mat magnitude, angle;
+    Mat hsv_split[3], hsv, rgb;
+    split(flow, flow_split);
+    cartToPolar(flow_split[0], flow_split[1], magnitude, angle, true);
+    normalize(magnitude, magnitude, 0, 1, NORM_MINMAX);
+    hsv_split[0] = angle; // already in degrees - no normalization needed
+    hsv_split[1] = Mat::ones(angle.size(), angle.type());
+    hsv_split[2] = magnitude;
+    merge(hsv_split, 3, hsv);
+    cvtColor(hsv, rgb, COLOR_HSV2BGR);
+    return rgb;
+}
+
+int main( int argc, char** argv )
+{
+    CommandLineParser parser(argc, argv, keys);
+    parser.about("OpenCV optical flow evaluation app");
+    if ( parser.has("help") || argc < 4 )
+    {
+        parser.printMessage();
+        printf("EXAMPLES:\n");
+        printf("./example_optflow_optical_flow_evaluation im1.png im2.png farneback -d \n");
+        printf("\t - compute flow field between im1 and im2 with farneback's method and display it");
+        printf("./example_optflow_optical_flow_evaluation im1.png im2.png simpleflow groundtruth.flo \n");
+        printf("\t - compute error statistics given the groundtruth; all pixels, endpoint error measure");
+        printf("./example_optflow_optical_flow_evaluation im1.png im2.png farneback groundtruth.flo -m=angular -r=untextured \n");
+        printf("\t - as before, but with changed error measure and stats computed only about \"untextured\" areas");
+        printf("\n\n Flow file format description: http://vision.middlebury.edu/flow/code/flow-code/README.txt\n\n");
+        return 0;
+    }
+    String i1_path = parser.get<String>(0);
+    String i2_path = parser.get<String>(1);
+    String method = parser.get<String>(2);
+    String groundtruth_path = parser.get<String>(3);
+    String error_measure = parser.get<String>("measure");
+    String region = parser.get<String>("region");
+    bool display_images = parser.has("display");
+
+    if ( !parser.check() )
+    {
+        parser.printErrors();
+        return 0;
+    }
+
+    Mat i1, i2;
+    Mat_<Point2f> flow, ground_truth;
+    Mat computed_errors;
+    i1 = imread(i1_path, 1);
+    i2 = imread(i2_path, 1);
+
+    if ( !i1.data || !i2.data )
+    {
+        printf("No image data \n");
+        return -1;
+    }
+    if ( i1.size() != i2.size() || i1.channels() != i2.channels() )
+    {
+        printf("Dimension mismatch between input images\n");
+        return -1;
+    }
+    // 8-bit images expected by all algorithms
+    if ( i1.depth() != CV_8U )
+        i1.convertTo(i1, CV_8U);
+    if ( i2.depth() != CV_8U )
+        i2.convertTo(i2, CV_8U);
+
+    if ( (method == "farneback" || method == "tvl1" || method == "deepflow") && i1.channels() == 3 )
+    {   // 1-channel images are expected
+        cvtColor(i1, i1, COLOR_BGR2GRAY);
+        cvtColor(i2, i2, COLOR_BGR2GRAY);
+    } else if ( method == "simpleflow" && i1.channels() == 1 )
+    {   // 3-channel images expected
+        cvtColor(i1, i1, COLOR_GRAY2BGR);
+        cvtColor(i2, i2, COLOR_GRAY2BGR);
+    }
+
+    flow = Mat(i1.size[0], i1.size[1], CV_32FC2);
+    Ptr<DenseOpticalFlow> algorithm;
+
+    if ( method == "farneback" )
+        algorithm = createOptFlow_Farneback();
+    else if ( method == "simpleflow" )
+        algorithm = createOptFlow_SimpleFlow();
+    else if ( method == "tvl1" )
+        algorithm = createOptFlow_DualTVL1();
+    else if ( method == "deepflow" )
+        algorithm = createOptFlow_DeepFlow();
+    else
+    {
+        printf("Wrong method!\n");
+        parser.printMessage();
+        return -1;
+    }
+
+    double startTick, time;
+    startTick = (double) getTickCount(); // measure time
+    algorithm->calc(i1, i2, flow);
+    time = ((double) getTickCount() - startTick) / getTickFrequency();
+    printf("\nTime [s]: %.3f\n", time);
+    if(display_images)
+    {
+        Mat flow_image = flowToDisplay(flow);
+        namedWindow( "Computed flow", WINDOW_AUTOSIZE );
+        imshow( "Computed flow", flow_image );
+    }
+
+    if ( !groundtruth_path.empty() )
+    { // compare to ground truth
+        ground_truth = optflow::readOpticalFlow(groundtruth_path);
+        if ( flow.size() != ground_truth.size() || flow.channels() != 2
+                || ground_truth.channels() != 2 )
+        {
+            printf("Dimension mismatch between the computed flow and the provided ground truth\n");
+            return -1;
+        }
+        if ( error_measure == "endpoint" )
+            computed_errors = endpointError(flow, ground_truth);
+        else if ( error_measure == "angular" )
+            computed_errors = angularError(flow, ground_truth);
+        else
+        {
+            printf("Invalid error measure! Available options: endpoint, angular\n");
+            return -1;
+        }
+
+        Mat mask;
+        if( region == "all" )
+            mask = Mat::ones(ground_truth.size(), CV_8U) * 255;
+        else if ( region == "discontinuities" )
+        {
+            Mat truth_merged, grad_x, grad_y, gradient;
+            vector<Mat> truth_split;
+            split(ground_truth, truth_split);
+            truth_merged = truth_split[0] + truth_split[1];
+
+            Sobel( truth_merged, grad_x, CV_16S, 1, 0, -1, 1, 0, BORDER_REPLICATE );
+            grad_x = abs(grad_x);
+            Sobel( truth_merged, grad_y, CV_16S, 0, 1, 1, 1, 0, BORDER_REPLICATE );
+            grad_y = abs(grad_y);
+            addWeighted(grad_x, 0.5, grad_y, 0.5, 0, gradient); //approximation!
+
+            Scalar s_mean;
+            s_mean = mean(gradient);
+            double threshold = s_mean[0]; // threshold value arbitrary
+            mask = gradient > threshold;
+            dilate(mask, mask, Mat::ones(9, 9, CV_8U));
+        }
+        else if ( region == "untextured" )
+        {
+            Mat i1_grayscale, grad_x, grad_y, gradient;
+            if( i1.channels() == 3 )
+                cvtColor(i1, i1_grayscale, COLOR_BGR2GRAY);
+            else
+                i1_grayscale = i1;
+            Sobel( i1_grayscale, grad_x, CV_16S, 1, 0, 7 );
+            grad_x = abs(grad_x);
+            Sobel( i1_grayscale, grad_y, CV_16S, 0, 1, 7 );
+            grad_y = abs(grad_y);
+            addWeighted(grad_x, 0.5, grad_y, 0.5, 0, gradient); //approximation!
+            GaussianBlur(gradient, gradient, Size(5,5), 1, 1);
+
+            Scalar s_mean;
+            s_mean = mean(gradient);
+            // arbitrary threshold value used - could be determined statistically from the image?
+            double threshold = 1000;
+            mask = gradient < threshold;
+            dilate(mask, mask, Mat::ones(3, 3, CV_8U));
+        }
+
+        else
+        {
+            printf("Invalid region selected! Available options: all, discontinuities, untextured");
+            return -1;
+        }
+
+        //masking out NaNs and incorrect GT values
+        Mat truth_split[2];
+        split(ground_truth, truth_split);
+        Mat abs_mask = Mat((abs(truth_split[0]) < 1e9) & (abs(truth_split[1]) < 1e9));
+        Mat nan_mask = Mat((truth_split[0]==truth_split[0]) & (truth_split[1] == truth_split[1]));
+        bitwise_and(abs_mask, nan_mask, nan_mask);
+
+        bitwise_and(nan_mask, mask, mask); //including the selected region
+
+        if(display_images) // display difference between computed and GT flow
+        {
+            Mat difference = ground_truth - flow;
+            Mat masked_difference;
+            difference.copyTo(masked_difference, mask);
+            Mat flow_image = flowToDisplay(masked_difference);
+            namedWindow( "Error map", WINDOW_AUTOSIZE );
+            imshow( "Error map", flow_image );
+        }
+
+        printf("Using %s error measure\n", error_measure.c_str());
+        calculateStats(computed_errors, mask, display_images);
+
+    }
+    if(display_images) // wait for the user to see all the images
+        waitKey(0);
+    return 0;
+
+}

modules/optflow/src/deepflow.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.
+ //
+ //
+ //                           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*/
+
+#include "precomp.hpp"
+#include <opencv2/highgui.hpp>
+
+namespace cv
+{
+namespace optflow
+{
+
+class OpticalFlowDeepFlow: public DenseOpticalFlow
+{
+public:
+    OpticalFlowDeepFlow();
+
+    void calc( InputArray I0, InputArray I1, InputOutputArray flow );
+    void collectGarbage();
+
+//    AlgorithmInfo* info() const;
+
+protected:
+    float sigma; // Gaussian smoothing parameter
+    int minSize; // minimal dimension of an image in the pyramid
+    float downscaleFactor; // scaling factor in the pyramid
+    int fixedPointIterations; // during each level of the pyramid
+    int sorIterations; // iterations of SOR
+    float alpha; // smoothness assumption weight
+    float delta; // color constancy weight
+    float gamma; // gradient constancy weight
+    float omega; // relaxation factor in SOR
+
+    float zeta; // added to the denomimnator of theta_0 (normaliation of the data term)
+    float epsilon; // robust penalizer const
+    int maxLayers; // max amount of layers in the pyramid
+
+private:
+    void calcOneLevel( const Mat I0, const Mat I1, Mat W );
+    Mat warpImage( const Mat input, const Mat flow );
+    void dataTerm( const Mat W, const Mat dW, const Mat Ix, const Mat Iy, const Mat Iz,
+            const Mat Ixx, const Mat Ixy, const Mat Iyy, const Mat Ixz, const Mat Iyz,
+            Mat a11, Mat a12, Mat a22, Mat b1, Mat b2 );
+    void smoothnessWeights( const Mat W, Mat weightsX, Mat weightsY );
+    void smoothnessTerm( const Mat W, const Mat weightsX, const Mat weightsY, Mat b1, Mat b2 );
+    void sorSolve( const Mat a11, const Mat a12, const Mat a22, const Mat b1, const Mat b2,
+            const Mat smoothX, const Mat smoothY, Mat dW );
+    void sorUnfolded( const Mat a11, const Mat a12, const Mat a22, const Mat b1, const Mat b2,
+            const Mat smoothX, const Mat smoothY, Mat dW );
+    std::vector<Mat> buildPyramid( const Mat& src );
+
+    int interpolationType;
+
+};
+
+OpticalFlowDeepFlow::OpticalFlowDeepFlow()
+{
+    // parameters
+    sigma = 0.6f;
+    minSize = 25;
+    downscaleFactor = 0.95f;
+    fixedPointIterations = 5;
+    sorIterations = 25;
+    alpha = 1.0f;
+    delta = 0.5f;
+    gamma = 5.0f;
+    omega = 1.6f;
+
+    //consts
+    interpolationType = INTER_LINEAR;
+    zeta = 0.1f;
+    epsilon = 0.001f;
+    maxLayers = 200;
+}
+
+std::vector<Mat> OpticalFlowDeepFlow::buildPyramid( const Mat& src )
+{
+    std::vector<Mat> pyramid;
+    pyramid.push_back(src);
+    Mat prev = pyramid[0];
+    int i = 0;
+    while ( i < this->maxLayers )
+    {
+        Mat next; //TODO: filtering at each level?
+        Size nextSize((int) (prev.cols * downscaleFactor + 0.5f),
+                        (int) (prev.rows * downscaleFactor + 0.5f));
+        if( nextSize.height <= minSize || nextSize.width <= minSize)
+            break;
+        resize(prev, next,
+                nextSize, 0, 0,
+                interpolationType);
+        pyramid.push_back(next);
+        prev = next;
+    }
+    return pyramid;
+}
+Mat OpticalFlowDeepFlow::warpImage( const Mat input, const Mat flow )
+{
+    // warps the image "backwards"
+    // if flow = computeFlow( I0, I1 ), then
+    // I0 = warpImage( I1, flow ) - approx.
+
+    Mat output;
+    Mat mapX = Mat(flow.size(), CV_32FC1);
+    Mat mapY = Mat(flow.size(), CV_32FC1);
+    const float *pFlow;
+    float *pMapX, *pMapY;
+    for ( int j = 0; j < flow.rows; ++j )
+    {
+        pFlow = flow.ptr<float>(j);
+        pMapX = mapX.ptr<float>(j);
+        pMapY = mapY.ptr<float>(j);
+        for ( int i = 0; i < flow.cols; ++i )
+        {
+            pMapX[i] = i + pFlow[2 * i];
+            pMapY[i] = j + pFlow[2 * i + 1];
+        }
+    }
+    remap(input, output, mapX, mapY, interpolationType, BORDER_TRANSPARENT);
+    return output;
+}
+void OpticalFlowDeepFlow::calc( InputArray _I0, InputArray _I1, InputOutputArray _flow )
+{
+    Mat I0temp = _I0.getMat();
+    Mat I1temp = _I1.getMat();
+
+    CV_Assert(I0temp.size() == I1temp.size());
+    CV_Assert(I0temp.type() == I1temp.type());
+    CV_Assert(I0temp.channels() == 1);
+    // TODO: currently only grayscale - data term could be computed in color version as well...
+
+    Mat I0, I1;
+
+    I0temp.convertTo(I0, CV_32F);
+    I1temp.convertTo(I1, CV_32F);
+
+    _flow.create(I0.size(), CV_32FC2);
+    Mat W = _flow.getMat(); // if any data present - will be discarded
+
+    // pre-smooth images
+    int kernelLen = ((int)floor(3 * sigma) * 2) + 1;
+    Size kernelSize(kernelLen, kernelLen);
+    GaussianBlur(I0, I0, kernelSize, sigma);
+    GaussianBlur(I1, I1, kernelSize, sigma);
+    // build down-sized pyramids
+    std::vector<Mat> pyramid_I0 = buildPyramid(I0);
+    std::vector<Mat> pyramid_I1 = buildPyramid(I1);
+    int levelCount = (int) pyramid_I0.size();
+
+    // initialize the first version of flow estimate to zeros
+    Size smallestSize = pyramid_I0[levelCount - 1].size();
+    W = Mat::zeros(smallestSize, CV_32FC2);
+
+    for ( int level = levelCount - 1; level >= 0; --level )
+    { //iterate through  all levels, beginning with the most coarse
+        calcOneLevel(pyramid_I0[level], pyramid_I1[level], W);
+        if ( level > 0 ) //not the last level
+        {
+            Mat temp;
+            Size newSize = pyramid_I0[level - 1].size();
+            resize(W, temp, newSize, 0, 0, interpolationType); //resize calculated flow
+            W = temp * (1.0f / downscaleFactor); //scale values
+        }
+    }
+    W.copyTo(_flow);
+}
+
+void OpticalFlowDeepFlow::calcOneLevel( const Mat I0, const Mat I1, Mat W )
+{
+    CV_DbgAssert( I0.size() == I1.size() );CV_DbgAssert( I0.type() == I1.type() );CV_DbgAssert( W.size() == I0.size() );
+
+    // linear equation systems
+    Size s = I0.size();
+    int t = CV_32F; // data type
+    Mat a11, a12, a22, b1, b2;
+    a11.create(s, t);
+    a12.create(s, t);
+    a22.create(s, t);
+    b1.create(s, t);
+    b2.create(s, t);
+    // diffusivity coeffs
+    Mat weightsX, weightsY;
+    weightsX.create(s, t);
+    weightsY.create(s, t);
+
+    Mat warpedI1 = warpImage(I1, W); // warped second image
+    Mat averageFrame = 0.5 * (I0 + warpedI1); // mean value of 2 frames - to compute derivatives on
+
+    //computing derivatives, notation as in Brox's paper
+    Mat Ix, Iy, Iz, Ixx, Ixy, Iyy, Ixz, Iyz;
+    int ddepth = -1; //as source image
+    int kernel_size = 1;
+
+    Sobel(averageFrame, Ix, ddepth, 1, 0, kernel_size, 1, 0.00, BORDER_REPLICATE);
+    Sobel(averageFrame, Iy, ddepth, 0, 1, kernel_size, 1, 0.00, BORDER_REPLICATE);
+    Iz.create(I1.size(), I1.type());
+    Iz = warpedI1 - I0;
+    Sobel(Ix, Ixx, ddepth, 1, 0, kernel_size, 1, 0.00, BORDER_REPLICATE);
+    Sobel(Ix, Ixy, ddepth, 0, 1, kernel_size, 1, 0.00, BORDER_REPLICATE);
+    Sobel(Iy, Iyy, ddepth, 0, 1, kernel_size, 1, 0.00, BORDER_REPLICATE);
+    Sobel(Iz, Ixz, ddepth, 1, 0, kernel_size, 1, 0.00, BORDER_REPLICATE);
+    Sobel(Iz, Iyz, ddepth, 0, 1, kernel_size, 1, 0.00, BORDER_REPLICATE);
+
+    Mat tempW = W.clone(); // flow version to be modified in each iteration
+    Mat dW = Mat::zeros(W.size(), W.type()); // flow increment
+
+    //fixed-point iterations
+    for ( int i = 0; i < fixedPointIterations; ++i )
+    {
+        dataTerm(W, dW, Ix, Iy, Iz, Ixx, Ixy, Iyy, Ixz, Iyz, a11, a12, a22, b1, b2);
+        smoothnessWeights(tempW, weightsX, weightsY);
+        smoothnessTerm(W, weightsX, weightsY, b1, b2);
+        sorSolve(a11, a12, a22, b1, b2, weightsX, weightsY, dW);
+        tempW = W + dW;
+    }
+    tempW.copyTo(W);
+}
+void OpticalFlowDeepFlow::dataTerm( const Mat W, const Mat dW, const Mat Ix, const Mat Iy,
+        const Mat Iz, const Mat Ixx, const Mat Ixy, const Mat Iyy, const Mat Ixz,
+        const Mat Iyz, Mat a11, Mat a12, Mat a22, Mat b1, Mat b2 )
+{
+    const float zeta_squared = zeta * zeta; // added in normalization factor to be non-zero
+    const float epsilon_squared = epsilon * epsilon;
+
+    const float *pIx, *pIy, *pIz;
+    const float *pIxx, *pIxy, *pIyy, *pIxz, *pIyz;
+    const float *pdU, *pdV; // accessing 2 layers of dW. Succesive columns interleave u and v
+    float *pa11, *pa12, *pa22, *pb1, *pb2; // linear equation sys. coeffs for each pixel
+
+    float derivNorm; //denominator of the spatial-derivative normalizing factor (theta_0)
+    float derivNorm2;
+    float Ik1z, Ik1zx, Ik1zy; // approximations of I^(k+1) values by Taylor expansions
+    float temp;
+    for ( int j = 0; j < W.rows; j++ ) //for each row
+    {
+        pIx = Ix.ptr<float>(j);
+        pIy = Iy.ptr<float>(j);
+        pIz = Iz.ptr<float>(j);
+        pIxx = Ixx.ptr<float>(j);
+        pIxy = Ixy.ptr<float>(j);
+        pIyy = Iyy.ptr<float>(j);
+        pIxz = Ixz.ptr<float>(j);
+        pIyz = Iyz.ptr<float>(j);
+
+        pa11 = a11.ptr<float>(j);
+        pa12 = a12.ptr<float>(j);
+        pa22 = a22.ptr<float>(j);
+        pb1 = b1.ptr<float>(j);
+        pb2 = b2.ptr<float>(j);
+
+        pdU = dW.ptr<float>(j);
+        pdV = pdU + 1;
+        for ( int i = 0; i < W.cols; i++ ) //for each pixel in the row
+        { // TODO: implement masking of points warped out of the image
+          //color constancy component
+            derivNorm = (*pIx) * (*pIx) + (*pIy) * (*pIy) + zeta_squared;
+            Ik1z = *pIz + (*pIx * *pdU) + (*pIy * *pdV);
+            temp = (0.5f*delta/3) / sqrt(Ik1z * Ik1z / derivNorm + epsilon_squared);
+            *pa11 = *pIx * *pIx * temp / derivNorm;
+            *pa12 = *pIx * *pIy * temp / derivNorm;
+            *pa22 = *pIy * *pIy * temp / derivNorm;
+            *pb1 = -*pIz * *pIx * temp / derivNorm;
+            *pb2 = -*pIz * *pIy * temp / derivNorm;
+
+            // gradient constancy component
+
+            derivNorm = *pIxx * *pIxx + *pIxy * *pIxy + zeta_squared;
+            derivNorm2 = *pIyy * *pIyy + *pIxy * *pIxy + zeta_squared;
+            Ik1zx = *pIxz + *pIxx * *pdU + *pIxy * *pdV;
+            Ik1zy = *pIyz + *pIxy * *pdU + *pIyy * *pdV;
+
+            temp = (0.5f*gamma/3)
+                    / sqrt(
+                            Ik1zx * Ik1zx / derivNorm + Ik1zy * Ik1zy / derivNorm2
+                                    + epsilon_squared);
+            *pa11 += temp * (*pIxx * *pIxx / derivNorm + *pIxy * *pIxy / derivNorm2);
+            *pa12 += temp * (*pIxx * *pIxy / derivNorm + *pIxy * *pIyy / derivNorm2);
+            *pa22 += temp * (*pIxy * *pIxy / derivNorm + *pIyy * *pIyy / derivNorm2);
+            *pb1 += -temp * (*pIxx * *pIxz / derivNorm + *pIxy * *pIyz / derivNorm2);
+            *pb2 += -temp * (*pIxy * *pIxz / derivNorm + *pIyy * *pIyz / derivNorm2);
+
+            ++pIx;
+            ++pIy;
+            ++pIz;
+            ++pIxx;
+            ++pIxy;
+            ++pIyy;
+            ++pIxz;
+            ++pIyz;
+            pdU += 2;
+            pdV += 2;
+            ++pa11;
+            ++pa12;
+            ++pa22;
+            ++pb1;
+            ++pb2;
+
+        }
+    }
+
+
+
+}
+void OpticalFlowDeepFlow::smoothnessWeights( const Mat W, Mat weightsX, Mat weightsY )
+{
+    float k[] = { -0.5, 0, 0.5 };
+    const float epsilon_squared = epsilon * epsilon;
+    Mat kernel_h = Mat(1, 3, CV_32FC1, k);
+    Mat kernel_v = Mat(3, 1, CV_32FC1, k);
+    Mat Wx, Wy; // partial derivatives of the flow
+    Mat S = Mat(W.size(), CV_32FC1); // sum of squared derivatives
+    weightsX = Mat::zeros(W.size(), CV_32FC1); //output - weights of smoothness terms in x and y directions
+    weightsY = Mat::zeros(W.size(), CV_32FC1);
+
+    filter2D(W, Wx, CV_32FC2, kernel_h);
+    filter2D(W, Wy, CV_32FC2, kernel_v);
+
+    const float * ux, *uy, *vx, *vy;
+    float * pS, *pWeight, *temp;
+
+    for ( int j = 0; j < S.rows; ++j )
+    {
+        ux = Wx.ptr<float>(j);
+        vx = ux + 1;
+        uy = Wy.ptr<float>(j);
+        vy = uy + 1;
+        pS = S.ptr<float>(j);
+        for ( int i = 0; i < S.cols; ++i )
+        {
+            *pS = alpha / sqrt(*ux * *ux + *vx * *vx + *uy * *uy + *vy * *vy + epsilon_squared);
+            ux += 2;
+            vx += 2;
+            uy += 2;
+            vy += 2;
+            ++pS;
+        }
+    }
+    // horizontal weights
+    for ( int j = 0; j < S.rows; ++j )
+    {
+        pWeight = weightsX.ptr<float>(j);
+        pS = S.ptr<float>(j);
+        for ( int i = 0; i < S.cols - 1; ++i )
+        {
+            *pWeight = *pS + *(pS + 1);
+            ++pS;
+            ++pWeight;
+        }
+    }
+    //vertical weights
+    for ( int j = 0; j < S.rows - 1; ++j )
+    {
+        pWeight = weightsY.ptr<float>(j);
+        pS = S.ptr<float>(j);
+        temp = S.ptr<float>(j + 1); // next row pointer for easy access
+        for ( int i = 0; i < S.cols; ++i )
+        {
+            *pWeight = *(pS++) + *(temp++);
+            ++pWeight;
+        }
+    }
+}
+void OpticalFlowDeepFlow::smoothnessTerm( const Mat W, const Mat weightsX, const Mat weightsY,
+        Mat b1, Mat b2 )
+{
+    float *pB1, *pB2;
+    const float *pU, *pV, *pWeight;
+    float iB1, iB2; // increments of b1 and b2
+    //horizontal direction - both U and V (b1 and b2)
+    for ( int j = 0; j < W.rows; j++ )
+    {
+        pB1 = b1.ptr<float>(j);
+        pB2 = b2.ptr<float>(j);
+        pU = W.ptr<float>(j);
+        pV = pU + 1;
+        pWeight = weightsX.ptr<float>(j);
+        for ( int i = 0; i < W.cols - 1; i++ )
+        {
+            iB1 = (*(pU + 2) - *pU) * *pWeight;
+            iB2 = (*(pV + 2) - *pV) * *pWeight;
+            *pB1 += iB1;
+            *(pB1 + 1) -= iB1;
+            *pB2 += iB2;
+            *(pB2 + 1) -= iB2;
+
+            pB1++;
+            pB2++;
+            pU += 2;
+            pV += 2;
+            pWeight++;
+        }
+    }
+    const float *pUnext, *pVnext; // temp pointers for next row
+    float *pB1next, *pB2next;
+    //vertical direction - both U and V
+    for ( int j = 0; j < W.rows - 1; j++ )
+    {
+        pB1 = b1.ptr<float>(j);
+        pB2 = b2.ptr<float>(j);
+        pU = W.ptr<float>(j);
+        pV = pU + 1;
+        pUnext = W.ptr<float>(j + 1);
+        pVnext = pUnext + 1;
+        pB1next = b1.ptr<float>(j + 1);
+        pB2next = b2.ptr<float>(j + 1);
+        pWeight = weightsY.ptr<float>(j);
+        for ( int i = 0; i < W.cols; i++ )
+        {
+            iB1 = (*pUnext - *pU) * *pWeight;
+            iB2 = (*pVnext - *pV) * *pWeight;
+            *pB1 += iB1;
+            *pB1next -= iB1;
+            *pB2 += iB2;
+            *pB2next -= iB2;
+
+            pB1++;
+            pB2++;
+            pU += 2;
+            pV += 2;
+            pWeight++;
+            pUnext += 2;
+            pVnext += 2;
+            pB1next++;
+            pB2next++;
+        }
+    }
+}
+
+void OpticalFlowDeepFlow::sorSolve( const Mat a11, const Mat a12, const Mat a22, const Mat b1,
+        const Mat b2, const Mat smoothX, const Mat smoothY, Mat dW )
+{
+    CV_Assert(a11.isContinuous());
+    CV_Assert(a12.isContinuous());
+    CV_Assert(a22.isContinuous());
+    CV_Assert(b1.isContinuous());
+    CV_Assert(b2.isContinuous());
+    CV_Assert(smoothX.isContinuous());
+    CV_Assert(smoothY.isContinuous());
+
+    if(dW.cols > 2 && dW.rows > 2)
+    {
+        sorUnfolded(a11, a12, a22, b1, b2, smoothX, smoothY, dW );
+        //more efficient version - this one is mostly for future reference and readability
+        return;
+    }
+    std::vector<Mat> dWChannels(2);
+    split(dW, dWChannels);
+
+    Mat *du = &(dWChannels[0]);
+    Mat *dv = &(dWChannels[1]);
+
+    CV_Assert(du->isContinuous());
+    CV_Assert(dv->isContinuous());
+
+    const float *pa11, *pa12, *pa22, *pb1, *pb2, *psmoothX, *psmoothY;
+    float *pdu, *pdv;
+    psmoothX = smoothX.ptr<float>(0);
+    psmoothY = smoothY.ptr<float>(0);
+    pdu = du->ptr<float>(0);
+    pdv = dv->ptr<float>(0);
+
+    float sigmaU, sigmaV, dPsi, A11, A22, A12, B1, B2, det;
+
+    int cols = dW.cols;
+    int rows = dW.rows;
+
+    int s = dW.cols; // step between rows
+
+    for ( int iter = 0; iter < sorIterations; ++iter )
+    {
+        pa11 = a11.ptr<float>(0);
+        pa12 = a12.ptr<float>(0);
+        pa22 = a22.ptr<float>(0);
+        pb1 = b1.ptr<float>(0);
+        pb2 = b2.ptr<float>(0);
+        for ( int j = 0; j < rows; ++j )
+        {
+            for ( int i = 0; i < cols; ++i )
+            {
+                int o = j * s + i;
+                if ( i == 0 && j == 0 )
+                {
+                    dPsi = psmoothX[o] + psmoothY[o];
+                    sigmaU = psmoothX[o] * pdu[o + 1] + psmoothY[o] * pdu[o + s];
+                    sigmaV = psmoothX[o] * pdv[o + 1] + psmoothY[o] * pdv[o + s];
+                } else if ( i == cols - 1 && j == 0 )
+                {
+                    dPsi = psmoothX[o - 1] + psmoothY[o];
+                    sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                            + psmoothY[o] * pdu[o + s];
+                    sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                            + psmoothY[o] * pdv[o + s];
+                } else if ( j == 0 )
+                {
+                    dPsi = psmoothX[o - 1] + psmoothX[o] + psmoothY[o];
+                    sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                            + psmoothX[o] * pdu[o + 1] + psmoothY[o] * pdu[o + s];
+                    sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                            + psmoothX[o] * pdv[o + 1] + psmoothY[o] * pdv[o + s];
+                } else if ( i == 0 && j == rows - 1 )
+                {
+                    dPsi = psmoothX[o] + psmoothY[o - s];
+                    sigmaU = psmoothX[o] * pdu[o + 1]
+                            + psmoothY[o - s] * pdu[o - s];
+                    sigmaV = psmoothX[o] * pdv[o + 1]
+                            + psmoothY[o - s] * pdv[o - s];
+                } else if ( i == cols - 1 && j == rows - 1 )
+                {
+                    dPsi = psmoothX[o - 1] + psmoothY[o - s];
+                    sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                            + psmoothY[o - s] * pdu[o - s];
+                    sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                            + psmoothY[o - s] * pdv[o - s];
+                } else if ( j == rows - 1 )
+                {
+                    dPsi = psmoothX[o - 1] + psmoothX[o] + psmoothY[o - s];
+                    sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                            + psmoothX[o] * pdu[o + 1]
+                            + psmoothY[o - s] * pdu[o - s];
+                    sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                            + psmoothX[o] * pdv[o + 1]
+                            + psmoothY[o - s] * pdv[o - s];
+                } else if ( i == 0 )
+                {
+                    dPsi = psmoothX[o] + psmoothY[o - s] + psmoothY[o];
+                    sigmaU = psmoothX[o] * pdu[o + 1]
+                            + psmoothY[o - s] * pdu[o - s]
+                            + psmoothY[o] * pdu[o + s];
+                    sigmaV = psmoothX[o] * pdv[o + 1]
+                            + psmoothY[o - s] * pdv[o - s]
+                            + psmoothY[o] * pdv[o + s];
+                } else if ( i == cols - 1 )
+                {
+                    dPsi = psmoothX[o - 1] + psmoothY[o - s] + psmoothY[o];
+                    sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                            + psmoothY[o - s] * pdu[o - s]
+                            + psmoothY[o] * pdu[o + s];
+                    sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                            + psmoothY[o - s] * pdv[o - s]
+                            + psmoothY[o] * pdv[o + s];
+                } else
+                {
+                    dPsi = psmoothX[o - 1] + psmoothX[o] + psmoothY[o - s]
+                            + psmoothY[o];
+                    sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                            + psmoothX[o] * pdu[o + 1]
+                            + psmoothY[o - s] * pdu[o - s]
+                            + psmoothY[o] * pdu[o + s];
+                    sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                            + psmoothX[o] * pdv[o + 1]
+                            + psmoothY[o - s] * pdv[o - s]
+                            + psmoothY[o] * pdv[o + s];
+                }
+                A11 = *pa22 + dPsi;
+                A12 = -*pa12;
+                A22 = *pa11 + dPsi;
+                det = A11 * A22 - A12 * A12;
+                A11 /= det;
+                A12 /= det;
+                A22 /= det;
+                B1 = *pb1 + sigmaU;
+                B2 = *pb2 + sigmaV;
+                pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+                pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+                ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+            }
+        }
+    }
+    merge(dWChannels, dW);
+}
+
+
+void OpticalFlowDeepFlow::sorUnfolded( const Mat a11, const Mat a12, const Mat a22, const Mat b1, const Mat b2,
+        const Mat smoothX, const Mat smoothY, Mat dW )
+{
+    // the same effect as sorSolve(), but written more efficiently
+    std::vector<Mat> dWChannels(2);
+    split(dW, dWChannels);
+
+    Mat *du = &(dWChannels[0]);
+    Mat *dv = &(dWChannels[1]);
+
+    CV_Assert(du->isContinuous());
+    CV_Assert(dv->isContinuous());
+
+    const float *pa11, *pa12, *pa22, *pb1, *pb2, *psmoothX, *psmoothY;
+    float *pdu, *pdv;
+
+
+    float sigmaU, sigmaV, dPsi, A11, A22, A12, B1, B2, det;
+
+    int cols = dW.cols;
+    int rows = dW.rows;
+
+    int s = dW.cols; // step between rows
+    int j, i, o; //row, column, offset
+
+    for ( int iter = 0; iter < sorIterations; ++iter )
+    {
+        pa11 = a11.ptr<float>(0);
+        pa12 = a12.ptr<float>(0);
+        pa22 = a22.ptr<float>(0);
+        pb1 = b1.ptr<float>(0);
+        pb2 = b2.ptr<float>(0);
+        psmoothX = smoothX.ptr<float>(0);
+        psmoothY = smoothY.ptr<float>(0);
+        pdu = du->ptr<float>(0);
+        pdv = dv->ptr<float>(0);
+
+        // first row
+        // first column
+        o=0;
+        dPsi = psmoothX[o] + psmoothY[o];
+        sigmaU = psmoothX[o] * pdu[o + 1] + psmoothY[o] * pdu[o + s];
+        sigmaV = psmoothX[o] * pdv[o + 1] + psmoothY[o] * pdv[o + s];
+        A11 = *pa22 + dPsi;
+        A12 = -*pa12;
+        A22 = *pa11 + dPsi;
+        det = A11 * A22 - A12 * A12;
+        A11 /= det;
+        A12 /= det;
+        A22 /= det;
+        B1 = *pb1 + sigmaU;
+        B2 = *pb2 + sigmaV;
+        pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+        pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+        ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+        // middle rows
+        for ( o = 1; o < cols-1; ++o )
+        {
+            dPsi = psmoothX[o - 1] + psmoothX[o] + psmoothY[o];
+            sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                    + psmoothX[o] * pdu[o + 1] + psmoothY[o] * pdu[o + s];
+            sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                    + psmoothX[o] * pdv[o + 1] + psmoothY[o] * pdv[o + s];
+            A11 = *pa22 + dPsi;
+            A12 = -*pa12;
+            A22 = *pa11 + dPsi;
+            det = A11 * A22 - A12 * A12;
+            A11 /= det;
+            A12 /= det;
+            A22 /= det;
+            B1 = *pb1 + sigmaU;
+            B2 = *pb2 + sigmaV;
+            pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+            pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+            ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+        }
+        // last column
+        dPsi = psmoothX[o - 1] + psmoothY[o];
+        sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                + psmoothY[o] * pdu[o + s];
+        sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                + psmoothY[o] * pdv[o + s];
+        A11 = *pa22 + dPsi;
+        A12 = -*pa12;
+        A22 = *pa11 + dPsi;
+        det = A11 * A22 - A12 * A12;
+        A11 /= det;
+        A12 /= det;
+        A22 /= det;
+        B1 = *pb1 + sigmaU;
+        B2 = *pb2 + sigmaV;
+        pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+        pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+        ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+        ++o;
+        //middle rows
+        for ( j = 1; j < rows - 1; ++j)
+        {
+            // first column
+            dPsi = psmoothX[o] + psmoothY[o - s] + psmoothY[o];
+            sigmaU = psmoothX[o] * pdu[o + 1]
+                    + psmoothY[o - s] * pdu[o - s]
+                    + psmoothY[o] * pdu[o + s];
+            sigmaV = psmoothX[o] * pdv[o + 1]
+                    + psmoothY[o - s] * pdv[o - s]
+                    + psmoothY[o] * pdv[o + s];
+            A11 = *pa22 + dPsi;
+            A12 = -*pa12;
+            A22 = *pa11 + dPsi;
+            det = A11 * A22 - A12 * A12;
+            A11 /= det;
+            A12 /= det;
+            A22 /= det;
+            B1 = *pb1 + sigmaU;
+            B2 = *pb2 + sigmaV;
+            pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+            pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+            ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+            ++o;
+            // middle columns
+            for ( i = 1; i < cols - 1; ++i)
+            {
+                dPsi = psmoothX[o - 1] + psmoothX[o] + psmoothY[o - s]
+                        + psmoothY[o];
+                sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                        + psmoothX[o] * pdu[o + 1]
+                        + psmoothY[o - s] * pdu[o - s]
+                        + psmoothY[o] * pdu[o + s];
+                sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                        + psmoothX[o] * pdv[o + 1]
+                        + psmoothY[o - s] * pdv[o - s]
+                        + psmoothY[o] * pdv[o + s];
+                A11 = *pa22 + dPsi;
+                A12 = -*pa12;
+                A22 = *pa11 + dPsi;
+                det = A11 * A22 - A12 * A12;
+                A11 /= det;
+                A12 /= det;
+                A22 /= det;
+                B1 = *pb1 + sigmaU;
+                B2 = *pb2 + sigmaV;
+                pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+                pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+                ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+                ++o;
+            }
+            //last column
+            dPsi = psmoothX[o - 1] + psmoothY[o - s] + psmoothY[o];
+            sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                    + psmoothY[o - s] * pdu[o - s]
+                    + psmoothY[o] * pdu[o + s];
+            sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                    + psmoothY[o - s] * pdv[o - s]
+                    + psmoothY[o] * pdv[o + s];
+            A11 = *pa22 + dPsi;
+            A12 = -*pa12;
+            A22 = *pa11 + dPsi;
+            det = A11 * A22 - A12 * A12;
+            A11 /= det;
+            A12 /= det;
+            A22 /= det;
+            B1 = *pb1 + sigmaU;
+            B2 = *pb2 + sigmaV;
+            pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+            pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+            ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+            ++o;
+        }
+        //last row
+        //first column
+        dPsi = psmoothX[o] + psmoothY[o - s];
+        sigmaU = psmoothX[o] * pdu[o + 1]
+                + psmoothY[o - s] * pdu[o - s];
+        sigmaV = psmoothX[o] * pdv[o + 1]
+                + psmoothY[o - s] * pdv[o - s];
+        A11 = *pa22 + dPsi;
+        A12 = -*pa12;
+        A22 = *pa11 + dPsi;
+        det = A11 * A22 - A12 * A12;
+        A11 /= det;
+        A12 /= det;
+        A22 /= det;
+        B1 = *pb1 + sigmaU;
+        B2 = *pb2 + sigmaV;
+        pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+        pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+        ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+        ++o;
+        //middle columns
+        for ( i = 1; i < cols - 1; ++i)
+        {
+            dPsi = psmoothX[o - 1] + psmoothX[o] + psmoothY[o - s];
+            sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                    + psmoothX[o] * pdu[o + 1]
+                    + psmoothY[o - s] * pdu[o - s];
+            sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                    + psmoothX[o] * pdv[o + 1]
+                    + psmoothY[o - s] * pdv[o - s];
+            A11 = *pa22 + dPsi;
+            A12 = -*pa12;
+            A22 = *pa11 + dPsi;
+            det = A11 * A22 - A12 * A12;
+            A11 /= det;
+            A12 /= det;
+            A22 /= det;
+            B1 = *pb1 + sigmaU;
+            B2 = *pb2 + sigmaV;
+            pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+            pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+            ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+            ++o;
+        }
+        //last column
+        dPsi = psmoothX[o - 1] + psmoothY[o - s];
+        sigmaU = psmoothX[o - 1] * pdu[o - 1]
+                + psmoothY[o - s] * pdu[o - s];
+        sigmaV = psmoothX[o - 1] * pdv[o - 1]
+                + psmoothY[o - s] * pdv[o - s];
+        A11 = *pa22 + dPsi;
+        A12 = -*pa12;
+        A22 = *pa11 + dPsi;
+        det = A11 * A22 - A12 * A12;
+        A11 /= det;
+        A12 /= det;
+        A22 /= det;
+        B1 = *pb1 + sigmaU;
+        B2 = *pb2 + sigmaV;
+        pdu[o] += omega * (A11 * B1 + A12 * B2 - pdu[o]);
+        pdv[o] += omega * (A12 * B1 + A22 * B2 - pdv[o]);
+        ++pa11; ++pa12; ++pa22; ++pb1; ++pb2;
+    }
+    merge(dWChannels, dW);
+
+}
+void OpticalFlowDeepFlow::collectGarbage()
+{
+
+}
+//
+//CV_INIT_ALGORITHM(OpticalFlowDeepFlow, "DenseOpticalFlow.DeepFlow",
+//        obj.info()->addParam(obj, "sigma", obj.sigma, false, 0, 0, "Gaussian blur parameter");
+//        obj.info()->addParam(obj, "alpha", obj.alpha, false, 0, 0, "Smoothness assumption weight");
+//        obj.info()->addParam(obj, "delta", obj.delta, false, 0, 0, "Color constancy weight");
+//        obj.info()->addParam(obj, "gamma", obj.gamma, false, 0, 0, "Gradient constancy weight");
+//        obj.info()->addParam(obj, "omega", obj.omega, false, 0, 0, "Relaxation factor in SOR");
+//        obj.info()->addParam(obj, "minSize", obj.minSize, false, 0, 0, "Min. image size in the pyramid");
+//        obj.info()->addParam(obj, "fixedPointIterations", obj.fixedPointIterations, false, 0, 0, "Fixed point iterations");
+//        obj.info()->addParam(obj, "sorIterations", obj.sorIterations, false, 0, 0, "SOR iterations");
+//        obj.info()->addParam(obj, "downscaleFactor", obj.downscaleFactor, false, 0, 0,"Downscale factor"))
+
+
+Ptr<DenseOpticalFlow> createOptFlow_DeepFlow()
+{
+    return makePtr<OpticalFlowDeepFlow>();
+}
+
+}//optflow
+}//cv

modules/optflow/src/interfaces.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.
+ //
+ //
+ //                           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*/
+#include "precomp.hpp"
+
+#include "opencv2/core.hpp"
+#include "opencv2/highgui.hpp"
+#include "opencv2/video.hpp"
+#include "opencv2/optflow.hpp"
+
+namespace cv
+{
+namespace optflow
+{
+class OpticalFlowSimpleFlow : public DenseOpticalFlow
+{
+public:
+    OpticalFlowSimpleFlow();
+    void calc(InputArray I0, InputArray I1, InputOutputArray flow);
+    void collectGarbage();
+//    AlgorithmInfo* info() const;
+
+protected:
+    int layers;
+    int averaging_radius;
+    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;
+};
+
+OpticalFlowSimpleFlow::OpticalFlowSimpleFlow()
+{
+    // values from the example app
+    layers = 3;
+    averaging_radius = 2;
+    max_flow = 4;
+    // values from the default function parameters
+    sigma_dist = 4.1;
+    sigma_color = 25.5;
+    postprocess_window = 18;
+    sigma_dist_fix = 55.0;
+    sigma_color_fix = 25.5;
+    occ_thr = 0.35;
+    upscale_averaging_radius = 18;
+    upscale_sigma_dist = 55.0;
+    upscale_sigma_color = 25.5;
+    speed_up_thr = 10;
+}
+
+void OpticalFlowSimpleFlow::calc(InputArray I0, InputArray I1, InputOutputArray flow)
+{
+    optflow::calcOpticalFlowSF(I0, I1, flow, layers, averaging_radius, max_flow, sigma_dist, sigma_color,
+                      postprocess_window, sigma_dist_fix, sigma_color_fix, occ_thr,
+                      upscale_averaging_radius, upscale_sigma_dist, upscale_sigma_color, speed_up_thr);
+}
+void OpticalFlowSimpleFlow::collectGarbage()
+{
+
+}
+//CV_INIT_ALGORITHM(OpticalFlowSimpleFlow, "DenseOpticalFlow.SimpleFlow"),
+//        obj.info()->addParam(obj, "layers", obj.layers);
+//        obj.info()->addParam(obj, "averaging_radius", obj.averaging_radius);
+//        obj.info()->addParam(obj, "max_flow", obj.max_flow);
+//        obj.info()->addParam(obj, "sigma_dist", obj.sigma_dist);
+//        obj.info()->addParam(obj, "sigma_color", obj.sigma_color);
+//        obj.info()->addParam(obj, "postprocess_window", obj.postprocess_window);
+//        obj.info()->addParam(obj, "sigma_dist_fix", obj.sigma_dist_fix);
+//        obj.info()->addParam(obj, "sigma_color_fix", obj.sigma_color_fix);
+//        obj.info()->addParam(obj, "occ_thr", obj.occ_thr);
+//        obj.info()->addParam(obj, "upscale_averaging_radius", obj.upscale_averaging_radius);
+//        obj.info()->addParam(obj, "upscale_sigma_dist", obj.upscale_sigma_dist);
+//        obj.info()->addParam(obj, "upscale_sigma_color", obj.upscale_sigma_color);
+//        obj.info()->addParam(obj, "speed_up_thr", obj.speed_up_thr))
+
+Ptr<DenseOpticalFlow> createOptFlow_SimpleFlow()
+{
+    return makePtr<OpticalFlowSimpleFlow>();
+}
+
+class OpticalFlowFarneback : public DenseOpticalFlow
+{
+public:
+    OpticalFlowFarneback();
+    void calc(InputArray I0, InputArray I1, InputOutputArray flow);
+    void collectGarbage();
+//    AlgorithmInfo* info() const;
+protected:
+    int numLevels;
+    double pyrScale;
+    bool fastPyramids;
+    int winSize;
+    int numIters;
+    int polyN;
+    double polySigma;
+    int flags;
+};
+
+OpticalFlowFarneback::OpticalFlowFarneback()
+{
+    // values copied from the FarnebackOpticalFlow class
+    numLevels = 5;
+    pyrScale = 0.5;
+    fastPyramids = false;
+    winSize = 13;
+    numIters = 10;
+    polyN = 5;
+    polySigma = 1.1;
+    flags = 0;
+}
+
+void OpticalFlowFarneback::calc(InputArray I0, InputArray I1, InputOutputArray flow)
+{
+    calcOpticalFlowFarneback(I0, I1, flow, pyrScale, numLevels, winSize, numIters, polyN, polySigma, flags);
+}
+void OpticalFlowFarneback::collectGarbage()
+{
+}
+
+//CV_INIT_ALGORITHM(OpticalFlowFarneback, "DenseOpticalFlow.Farneback",
+//        obj.info()->addParam(obj, "numLevels", obj.numLevels);
+//        obj.info()->addParam(obj, "pyrScale", obj.pyrScale);
+//        obj.info()->addParam(obj, "fastPyramids", obj.fastPyramids);
+//        obj.info()->addParam(obj, "winSize", obj.winSize);
+//        obj.info()->addParam(obj, "numIters", obj.numIters);
+//        obj.info()->addParam(obj, "polyN", obj.polyN);
+//        obj.info()->addParam(obj, "polySigma", obj.polySigma);
+//        obj.info()->addParam(obj, "flags", obj.flags))
+
+
+
+
+Ptr<DenseOpticalFlow> createOptFlow_Farneback()
+{
+    return makePtr<OpticalFlowFarneback>();
+}
+}
+}

modules/optflow/src/optical_flow_io.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.
+ //
+ //
+ //                           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*/
+#include "precomp.hpp"
+#include<iostream>
+#include<fstream>
+
+namespace cv {
+namespace optflow {
+const float FLOW_TAG_FLOAT = 202021.25f;
+const char *FLOW_TAG_STRING = "PIEH";
+CV_EXPORTS_W Mat readOpticalFlow( const String& path )
+{
+//    CV_Assert(sizeof(float) == 4);
+    //FIXME: ensure right sizes of int and float - here and in writeOpticalFlow()
+
+    Mat_<Point2f> flow;
+    std::ifstream file(path.c_str(), std::ios_base::binary);
+    if ( !file.good() )
+        return flow; // no file - return empty matrix
+
+    float tag;
+    file.read((char*) &tag, sizeof(float));
+    if ( tag != FLOW_TAG_FLOAT )
+        return flow;
+
+    int width, height;
+
+    file.read((char*) &width, 4);
+    file.read((char*) &height, 4);
+
+    flow.create(height, width);
+
+    for ( int i = 0; i < flow.rows; ++i )
+    {
+        for ( int j = 0; j < flow.cols; ++j )
+        {
+            Point2f u;
+            file.read((char*) &u.x, sizeof(float));
+            file.read((char*) &u.y, sizeof(float));
+            if ( !file.good() )
+            {
+                flow.release();
+                return flow;
+            }
+
+            flow(i, j) = u;
+        }
+    }
+    file.close();
+    return flow;
+}
+CV_EXPORTS_W bool writeOpticalFlow( const String& path, InputArray flow )
+{
+//    CV_Assert(sizeof(float) == 4);
+
+    const int nChannels = 2;
+
+    Mat input = flow.getMat();
+    if ( input.channels() != nChannels || input.depth() != CV_32F || path.length() == 0 )
+        return false;
+
+    std::ofstream file(path.c_str(), std::ofstream::binary);
+    if ( !file.good() )
+        return false;
+
+    int nRows, nCols;
+
+    nRows = (int) input.size().height;
+    nCols = (int) input.size().width;
+
+    const int headerSize = 12;
+    char header[headerSize];
+    memcpy(header, FLOW_TAG_STRING, 4);
+    // size of ints is known - has been asserted in the current function
+    memcpy(header + 4, reinterpret_cast<const char*>(&nCols), sizeof(nCols));
+    memcpy(header + 8, reinterpret_cast<const char*>(&nRows), sizeof(nRows));
+    file.write(header, headerSize);
+    if ( !file.good() )
+        return false;
+
+//    if ( input.isContinuous() ) //matrix is continous - treat it as a single row
+//    {
+//        nCols *= nRows;
+//        nRows = 1;
+//    }
+
+    int row;
+    char* p;
+    for ( row = 0; row < nRows; row++ )
+    {
+        p = input.ptr<char>(row);
+        file.write(p, nCols * nChannels * sizeof(float));
+        if ( !file.good() )
+            return false;
+    }
+    file.close();
+    return true;
+}
+}
+}