Moved deepflow to contrib, minor improvements

modules/optflow/include/opencv2/optflow.hpp

@@ -70,7 +70,7 @@ 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();
+CV_EXPORTS_W Ptr<DenseOpticalFlow> createOptFlow_DeepFlow();
 
 // Additional interface to the SimpleFlow algorithm - calcOpticalFlowSF()
 CV_EXPORTS_W Ptr<DenseOpticalFlow> createOptFlow_SimpleFlow();

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 );
+    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 = (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 = 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.5*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.5*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());
+
+    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];
+                }
+                A22 = *pa11 + dPsi;
+                A11 = *pa22 + dPsi;
+                A12 = -*pa12;
+                det = A11 * A22 - A12 * A12;
+                A22 /= det;
+                A11 /= det;
+                A12 /= 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