diff --git a/modules/optflow/doc/dense_optflow.rst b/modules/optflow/doc/dense_optflow.rst index 17c1eb99d..0a678a7fc 100644 --- a/modules/optflow/doc/dense_optflow.rst +++ b/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 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. diff --git a/modules/optflow/doc/optflow.rst b/modules/optflow/doc/optflow.rst index 56f3366df..5b637167e 100644 --- a/modules/optflow/doc/optflow.rst +++ b/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 diff --git a/modules/optflow/doc/optflow_io.rst b/modules/optflow/doc/optflow_io.rst new file mode 100644 index 000000000..ed3319970 --- /dev/null +++ b/modules/optflow/doc/optflow_io.rst @@ -0,0 +1,36 @@ +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 diff --git a/modules/optflow/include/opencv2/optflow.hpp b/modules/optflow/include/opencv2/optflow.hpp index 7f62f7406..e7427d95a 100644 --- a/modules/optflow/include/opencv2/optflow.hpp +++ b/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 createOptFlow_DeepFlow(); + +// Additional interface to the SimpleFlow algorithm - calcOpticalFlowSF() +CV_EXPORTS_W Ptr createOptFlow_SimpleFlow(); + +// Additional interface to the Farneback's algorithm - calcOpticalFlowFarneback() +CV_EXPORTS_W Ptr createOptFlow_Farneback(); +} //optflow } #include "opencv2/optflow/motempl.hpp" diff --git a/modules/optflow/samples/optical_flow_evaluation.cpp b/modules/optflow/samples/optical_flow_evaluation.cpp new file mode 100644 index 000000000..c65412a3e --- /dev/null +++ b/modules/optflow/samples/optical_flow_evaluation.cpp @@ -0,0 +1,367 @@ +#include "opencv2/highgui.hpp" +#include "opencv2/video.hpp" +#include "opencv2/optflow.hpp" +#include + +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_& flow1, const Mat_& 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(i, j) = sqrt((float)diff.ddot(diff)); //distance + } else + result.at(i, j) = NAN; + } + } + return result; +} +static Mat angularError( const Mat_& flow1, const Mat_& 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(i, j) = acos((float)(u1.ddot(u2) / norm(u1) * norm(u2))); + else + result.at(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(i, j) != 0 ) + { + ++all; + if ( errors.at(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(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(0); + String i2_path = parser.get(1); + String method = parser.get(2); + String groundtruth_path = parser.get(3); + String error_measure = parser.get("measure"); + String region = parser.get("region"); + bool display_images = parser.has("display"); + + if ( !parser.check() ) + { + parser.printErrors(); + return 0; + } + + Mat i1, i2; + Mat_ 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 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 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; + +} diff --git a/modules/optflow/src/deepflow.cpp b/modules/optflow/src/deepflow.cpp new file mode 100644 index 000000000..334b1d2b9 --- /dev/null +++ b/modules/optflow/src/deepflow.cpp @@ -0,0 +1,868 @@ +/*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 + +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 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 OpticalFlowDeepFlow::buildPyramid( const Mat& src ) +{ + std::vector 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(j); + pMapX = mapX.ptr(j); + pMapY = mapY.ptr(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 pyramid_I0 = buildPyramid(I0); + std::vector 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(j); + pIy = Iy.ptr(j); + pIz = Iz.ptr(j); + pIxx = Ixx.ptr(j); + pIxy = Ixy.ptr(j); + pIyy = Iyy.ptr(j); + pIxz = Ixz.ptr(j); + pIyz = Iyz.ptr(j); + + pa11 = a11.ptr(j); + pa12 = a12.ptr(j); + pa22 = a22.ptr(j); + pb1 = b1.ptr(j); + pb2 = b2.ptr(j); + + pdU = dW.ptr(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(j); + vx = ux + 1; + uy = Wy.ptr(j); + vy = uy + 1; + pS = S.ptr(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(j); + pS = S.ptr(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(j); + pS = S.ptr(j); + temp = S.ptr(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(j); + pB2 = b2.ptr(j); + pU = W.ptr(j); + pV = pU + 1; + pWeight = weightsX.ptr(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(j); + pB2 = b2.ptr(j); + pU = W.ptr(j); + pV = pU + 1; + pUnext = W.ptr(j + 1); + pVnext = pUnext + 1; + pB1next = b1.ptr(j + 1); + pB2next = b2.ptr(j + 1); + pWeight = weightsY.ptr(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 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(0); + psmoothY = smoothY.ptr(0); + pdu = du->ptr(0); + pdv = dv->ptr(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(0); + pa12 = a12.ptr(0); + pa22 = a22.ptr(0); + pb1 = b1.ptr(0); + pb2 = b2.ptr(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 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(0); + pa12 = a12.ptr(0); + pa22 = a22.ptr(0); + pb1 = b1.ptr(0); + pb2 = b2.ptr(0); + psmoothX = smoothX.ptr(0); + psmoothY = smoothY.ptr(0); + pdu = du->ptr(0); + pdv = dv->ptr(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 createOptFlow_DeepFlow() +{ + return makePtr(); +} + +}//optflow +}//cv diff --git a/modules/optflow/src/interfaces.cpp b/modules/optflow/src/interfaces.cpp new file mode 100644 index 000000000..f11f54ed0 --- /dev/null +++ b/modules/optflow/src/interfaces.cpp @@ -0,0 +1,183 @@ +/*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 createOptFlow_SimpleFlow() +{ + return makePtr(); +} + +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 createOptFlow_Farneback() +{ + return makePtr(); +} +} +} diff --git a/modules/optflow/src/optical_flow_io.cpp b/modules/optflow/src/optical_flow_io.cpp new file mode 100644 index 000000000..ad70b121f --- /dev/null +++ b/modules/optflow/src/optical_flow_io.cpp @@ -0,0 +1,139 @@ +/*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 +#include + +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_ 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(&nCols), sizeof(nCols)); + memcpy(header + 8, reinterpret_cast(&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(row); + file.write(p, nCols * nChannels * sizeof(float)); + if ( !file.good() ) + return false; + } + file.close(); + return true; +} +} +}