Merge pull request #710 from VladX:optflow

modules/optflow/include/opencv2/optflow.hpp

@@ -43,6 +43,9 @@ the use of this software, even if advised of the possibility of such damage.
 #include "opencv2/core.hpp"
 #include "opencv2/video.hpp"
 
+#include "opencv2/optflow/pcaflow.hpp"
+#include "opencv2/optflow/sparse_matching_gpc.hpp"
+
 /**
 @defgroup optflow Optical Flow Algorithms
 

modules/optflow/include/opencv2/optflow/pcaflow.hpp

@@ -0,0 +1,134 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2016, OpenCV Foundation, 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:
+
+  * Redistributions of source code must retain the above copyright notice,
+    this list of conditions and the following disclaimer.
+
+  * Redistributions 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.
+
+  * Neither the names of the copyright holders nor the names of the contributors
+    may 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 copyright holders 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.
+*/
+
+/*
+Implementation of the PCAFlow algorithm from the following paper:
+http://files.is.tue.mpg.de/black/papers/cvpr2015_pcaflow.pdf
+
+@inproceedings{Wulff:CVPR:2015,
+  title = {Efficient Sparse-to-Dense Optical Flow Estimation using a Learned Basis and Layers},
+  author = {Wulff, Jonas and Black, Michael J.},
+  booktitle = { IEEE Conf. on Computer Vision and Pattern Recognition (CVPR) 2015},
+  month = jun,
+  year = {2015}
+}
+
+There are some key differences which distinguish this algorithm from the original PCAFlow (see paper):
+  - Discrete Cosine Transform basis is used instead of basis extracted with PCA.
+    Reasoning: DCT basis has comparable performance and it doesn't require additional storage space.
+    Also, this decision helps to avoid overloading the algorithm with a lot of external input.
+  - Usage of built-in OpenCV feature tracking instead of libviso.
+*/
+
+#ifndef __OPENCV_OPTFLOW_PCAFLOW_HPP__
+#define __OPENCV_OPTFLOW_PCAFLOW_HPP__
+
+#include "opencv2/core.hpp"
+#include "opencv2/video.hpp"
+
+namespace cv
+{
+namespace optflow
+{
+
+/*
+ * This class can be used for imposing a learned prior on the resulting optical flow.
+ * Solution will be regularized according to this prior.
+ * You need to generate appropriate prior file with "learn_prior.py" script beforehand.
+ */
+class CV_EXPORTS_W PCAPrior
+{
+private:
+  Mat L1;
+  Mat L2;
+  Mat c1;
+  Mat c2;
+
+public:
+  PCAPrior( const char *pathToPrior );
+
+  int getPadding() const { return L1.size().height; }
+
+  int getBasisSize() const { return L1.size().width; }
+
+  void fillConstraints( float *A1, float *A2, float *b1, float *b2 ) const;
+};
+
+class CV_EXPORTS_W OpticalFlowPCAFlow : public DenseOpticalFlow
+{
+protected:
+  const Ptr<const PCAPrior> prior;
+  const Size basisSize;
+  const float sparseRate;              // (0 .. 0.1)
+  const float retainedCornersFraction; // [0 .. 1]
+  const float occlusionsThreshold;
+  const float dampingFactor;
+  const float claheClip;
+  bool useOpenCL;
+
+public:
+  OpticalFlowPCAFlow( Ptr<const PCAPrior> _prior = Ptr<const PCAPrior>(), const Size _basisSize = Size( 18, 14 ),
+                      float _sparseRate = 0.024, float _retainedCornersFraction = 0.2,
+                      float _occlusionsThreshold = 0.0003, float _dampingFactor = 0.00002, float _claheClip = 14 );
+
+  void calc( InputArray I0, InputArray I1, InputOutputArray flow );
+  void collectGarbage();
+
+private:
+  void findSparseFeatures( UMat &from, UMat &to, std::vector<Point2f> &features,
+                           std::vector<Point2f> &predictedFeatures ) const;
+
+  void removeOcclusions( UMat &from, UMat &to, std::vector<Point2f> &features,
+                         std::vector<Point2f> &predictedFeatures ) const;
+
+  void getSystem( OutputArray AOut, OutputArray b1Out, OutputArray b2Out, const std::vector<Point2f> &features,
+                  const std::vector<Point2f> &predictedFeatures, const Size size );
+
+  void getSystem( OutputArray A1Out, OutputArray A2Out, OutputArray b1Out, OutputArray b2Out,
+                  const std::vector<Point2f> &features, const std::vector<Point2f> &predictedFeatures,
+                  const Size size );
+
+  OpticalFlowPCAFlow& operator=( const OpticalFlowPCAFlow& ); // make it non-assignable
+};
+
+CV_EXPORTS_W Ptr<DenseOpticalFlow> createOptFlow_PCAFlow();
+}
+}
+
+#endif

modules/optflow/include/opencv2/optflow/sparse_matching_gpc.hpp

@@ -0,0 +1,184 @@
+/*
+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
+                       (3-clause BSD License)
+
+Copyright (C) 2016, OpenCV Foundation, 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:
+
+  * Redistributions of source code must retain the above copyright notice,
+    this list of conditions and the following disclaimer.
+
+  * Redistributions 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.
+
+  * Neither the names of the copyright holders nor the names of the contributors
+    may 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 copyright holders 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.
+*/
+
+/*
+Implementation of the Global Patch Collider algorithm from the following paper:
+http://research.microsoft.com/en-us/um/people/pkohli/papers/wfrik_cvpr2016.pdf
+
+@InProceedings{Wang_2016_CVPR,
+ author = {Wang, Shenlong and Ryan Fanello, Sean and Rhemann, Christoph and Izadi, Shahram and Kohli, Pushmeet},
+ title = {The Global Patch Collider},
+ booktitle = {The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)},
+ month = {June},
+ year = {2016}
+}
+*/
+
+#ifndef __OPENCV_OPTFLOW_SPARSE_MATCHING_GPC_HPP__
+#define __OPENCV_OPTFLOW_SPARSE_MATCHING_GPC_HPP__
+
+#include "opencv2/core.hpp"
+
+namespace cv
+{
+namespace optflow
+{
+
+struct CV_EXPORTS_W GPCPatchDescriptor
+{
+  static const unsigned nFeatures = 18; // number of features in a patch descriptor
+  Vec< double, nFeatures > feature;
+
+  GPCPatchDescriptor( const Mat *imgCh, int i, int j );
+};
+
+typedef std::pair< GPCPatchDescriptor, GPCPatchDescriptor > GPCPatchSample;
+typedef std::vector< GPCPatchSample > GPCSamplesVector;
+
+/** @brief Class encapsulating training samples.
+ */
+class CV_EXPORTS_W GPCTrainingSamples
+{
+private:
+  GPCSamplesVector samples;
+
+public:
+  /** @brief This function can be used to extract samples from a pair of images and a ground truth flow.
+   * Sizes of all the provided vectors must be equal.
+   */
+  static Ptr< GPCTrainingSamples > create( const std::vector< String > &imagesFrom, const std::vector< String > &imagesTo,
+                                           const std::vector< String > &gt );
+
+  size_t size() const { return samples.size(); }
+
+  operator GPCSamplesVector() const { return samples; }
+
+  operator GPCSamplesVector &() { return samples; }
+};
+
+class CV_EXPORTS_W GPCTree : public Algorithm
+{
+public:
+  struct Node
+  {
+    Vec< double, GPCPatchDescriptor::nFeatures > coef; // hyperplane coefficients
+    double rhs;
+    unsigned left;
+    unsigned right;
+
+    bool operator==( const Node &n ) const { return coef == n.coef && rhs == n.rhs && left == n.left && right == n.right; }
+  };
+
+private:
+  typedef GPCSamplesVector::iterator SIter;
+
+  std::vector< Node > nodes;
+
+  bool trainNode( size_t nodeId, SIter begin, SIter end, unsigned depth );
+
+public:
+  void train( GPCSamplesVector &samples );
+
+  void write( FileStorage &fs ) const;
+
+  void read( const FileNode &fn );
+
+  static Ptr< GPCTree > create() { return makePtr< GPCTree >(); }
+
+  bool operator==( const GPCTree &t ) const { return nodes == t.nodes; }
+};
+
+template < int T > class CV_EXPORTS_W GPCForest : public Algorithm
+{
+private:
+  GPCTree tree[T];
+
+public:
+  /** @brief Train the forest using one sample set for every tree.
+   * Please, consider using the next method instead of this one for better quality.
+   */
+  void train( GPCSamplesVector &samples )
+  {
+    for ( int i = 0; i < T; ++i )
+      tree[i].train( samples );
+  }
+
+  /** @brief Train the forest using individual samples for each tree.
+   * It is generally better to use this instead of the first method.
+   */
+  void train( const std::vector< String > &imagesFrom, const std::vector< String > &imagesTo, const std::vector< String > &gt )
+  {
+    for ( int i = 0; i < T; ++i )
+    {
+      Ptr< GPCTrainingSamples > samples = GPCTrainingSamples::create( imagesFrom, imagesTo, gt ); // Create training set for the tree
+      tree[i].train( *samples );
+    }
+  }
+
+  void write( FileStorage &fs ) const
+  {
+    fs << "ntrees" << T << "trees"
+       << "[";
+    for ( int i = 0; i < T; ++i )
+    {
+      fs << "{";
+      tree[i].write( fs );
+      fs << "}";
+    }
+    fs << "]";
+  }
+
+  void read( const FileNode &fn )
+  {
+    CV_Assert( T == (int)fn["ntrees"] );
+    FileNodeIterator it = fn["trees"].begin();
+    for ( int i = 0; i < T; ++i, ++it )
+      tree[i].read( *it );
+  }
+
+  static Ptr< GPCForest > create() { return makePtr< GPCForest >(); }
+};
+}
+
+CV_EXPORTS void write( FileStorage &fs, const String &name, const optflow::GPCTree::Node &node );
+
+CV_EXPORTS void read( const FileNode &fn, optflow::GPCTree::Node &node, optflow::GPCTree::Node );
+}
+
+#endif

modules/optflow/samples/gpc_train.cpp

@@ -0,0 +1,31 @@
+#include "opencv2/optflow.hpp"
+#include <iostream>
+
+const int nTrees = 5;
+
+int main( int argc, const char **argv )
+{
+  int nSequences = argc - 1;
+
+  if ( nSequences <= 0 || nSequences % 3 != 0 )
+  {
+    std::cerr << "Usage: " << argv[0] << " ImageFrom1 ImageTo1 GroundTruth1 ... ImageFromN ImageToN GroundTruthN" << std::endl;
+    return 1;
+  }
+
+  nSequences /= 3;
+  std::vector< cv::String > img1, img2, gt;
+
+  for ( int i = 0; i < nSequences; ++i )
+  {
+    img1.push_back( argv[1 + i * 3] );
+    img2.push_back( argv[1 + i * 3 + 1] );
+    gt.push_back( argv[1 + i * 3 + 2] );
+  }
+
+  cv::Ptr< cv::optflow::GPCForest< nTrees > > forest = cv::optflow::GPCForest< nTrees >::create();
+  forest->train( img1, img2, gt );
+  forest->save( "forest.dump" );
+
+  return 0;
+}

modules/optflow/samples/optical_flow_evaluation.cpp

@@ -1,6 +1,7 @@
 #include "opencv2/highgui.hpp"
 #include "opencv2/video.hpp"
 #include "opencv2/optflow.hpp"
+#include "opencv2/core/ocl.hpp"
 #include <fstream>
 #include <limits>
 
@@ -11,11 +12,13 @@ using namespace optflow;
 const String keys = "{help h usage ? |      | print this message   }"
         "{@image1        |      | image1               }"
         "{@image2        |      | image2               }"
-        "{@algorithm     |      | [farneback, simpleflow, tvl1, deepflow, sparsetodenseflow, DISflow_ultrafast, DISflow_fast, DISflow_medium] }"
+        "{@algorithm     |      | [farneback, simpleflow, tvl1, deepflow, sparsetodenseflow, pcaflow, DISflow_ultrafast, DISflow_fast, DISflow_medium] }"
         "{@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) }";
+        "{d display      |      | display additional info images (pauses program execution) }"
+        "{g gpu          |      | use OpenCL}"
+        "{prior          |      | path to a prior file for PCAFlow}";
 
 inline bool isFlowCorrect( const Point2f u )
 {
@@ -200,6 +203,7 @@ int main( int argc, char** argv )
     String error_measure = parser.get<String>("measure");
     String region = parser.get<String>("region");
     bool display_images = parser.has("display");
+    const bool useGpu = parser.has("gpu");
 
     if ( !parser.check() )
     {
@@ -207,6 +211,9 @@ int main( int argc, char** argv )
         return 0;
     }
 
+    cv::ocl::setUseOpenCL(useGpu);
+    printf("OpenCL Enabled: %u\n", useGpu && cv::ocl::haveOpenCL());
+
     Mat i1, i2;
     Mat_<Point2f> flow, ground_truth;
     Mat computed_errors;
@@ -252,6 +259,15 @@ int main( int argc, char** argv )
         algorithm = createOptFlow_DeepFlow();
     else if ( method == "sparsetodenseflow" )
         algorithm = createOptFlow_SparseToDense();
+    else if ( method == "pcaflow" ) {
+        if ( parser.has("prior") ) {
+            String prior = parser.get<String>("prior");
+            printf("Using prior file: %s\n", prior.c_str());
+            algorithm = makePtr<OpticalFlowPCAFlow>(makePtr<PCAPrior>(prior.c_str()));
+        }
+        else
+            algorithm = createOptFlow_PCAFlow();
+    }
     else if ( method == "DISflow_ultrafast" )
         algorithm = createOptFlow_DIS(DISOpticalFlow::PRESET_ULTRAFAST);
     else if (method == "DISflow_fast")
@@ -267,7 +283,12 @@ int main( int argc, char** argv )
 
     double startTick, time;
     startTick = (double) getTickCount(); // measure time
-    algorithm->calc(i1, i2, flow);
+
+    if (useGpu)
+        algorithm->calc(i1, i2, flow.getUMat(ACCESS_RW));
+    else
+        algorithm->calc(i1, i2, flow);
+
     time = ((double) getTickCount() - startTick) / getTickFrequency();
     printf("\nTime [s]: %.3f\n", time);
     if(display_images)

modules/optflow/src/learn_prior.py

@@ -0,0 +1,166 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import numpy as np
+import cv2
+import struct
+import argparse
+from math import sqrt
+
+argparser = argparse.ArgumentParser(
+    description='''Use this script to generate prior for using with PCAFlow.
+Basis size here must match corresponding parameter in the PCAFlow.
+Gamma should be selected experimentally.''')
+
+argparser.add_argument('-f',
+                       '--files',
+                       nargs='+',
+                       help='List of optical flow .flo files for learning. You can pass a directory here and it will be scanned recursively for .flo files.',
+                       required=True)
+argparser.add_argument('-o',
+                       '--output',
+                       help='Output file for prior',
+                       required=True)
+argparser.add_argument('--width',
+                       type=int,
+                       help='Size of the basis first dimension',
+                       required=True,
+                       default=18)
+argparser.add_argument('--height',
+                       type=int,
+                       help='Size of the basis second dimension',
+                       required=True,
+                       default=14)
+argparser.add_argument(
+    '-g',
+    '--gamma',
+    type=float,
+    help='Amount of regularization. The greater this parameter, the bigger will be an impact of the regularization.',
+    required=True)
+args = argparser.parse_args()
+
+basis_size = (args.height, args.width)
+gamma = args.gamma
+
+
+def find_flo(pp):
+    f = []
+    for p in pp:
+        if os.path.isfile(p):
+            f.append(p)
+        else:
+            for root, subdirs, files in os.walk(p):
+                f += map(lambda x: os.path.join(root, x),
+                         filter(lambda x: x.split('.')[-1] == 'flo', files))
+    return list(set(f))
+
+
+def load_flo(flo):
+    with open(flo, 'rb') as f:
+        magic = np.fromfile(f, np.float32, count=1)[0]
+        if 202021.25 != magic:
+            print('Magic number incorrect. Invalid .flo file')
+        else:
+            w = np.fromfile(f, np.int32, count=1)[0]
+            h = np.fromfile(f, np.int32, count=1)[0]
+            print('Reading %dx%d flo file %s' % (w, h, flo))
+            data = np.fromfile(f, np.float32, count=2 * w * h)
+            # Reshape data into 3D array (columns, rows, bands)
+            flow = np.reshape(data, (h, w, 2))
+            return flow[:, :, 0], flow[:, :, 1]
+
+
+def get_w(m):
+    s = m.shape
+    w = cv2.dct(m)
+    w *= 2.0 / sqrt(s[0] * s[1])
+    #w[0,0] *= 0.5
+    w[:, 0] *= sqrt(0.5)
+    w[0, :] *= sqrt(0.5)
+    w = w[0:basis_size[0], 0:basis_size[1]].transpose().flatten()
+    return w
+
+
+w1 = []
+w2 = []
+
+for flo in find_flo(args.files):
+    x, y = load_flo(flo)
+    w1.append(get_w(x))
+    w2.append(get_w(y))
+
+w1mean = sum(w1) / len(w1)
+w2mean = sum(w2) / len(w2)
+
+for i in xrange(len(w1)):
+    w1[i] -= w1mean
+for i in xrange(len(w2)):
+    w2[i] -= w2mean
+
+Q1 = sum([w1[i].reshape(-1, 1).dot(w1[i].reshape(1, -1))
+          for i in xrange(len(w1))]) / len(w1)
+Q2 = sum([w2[i].reshape(-1, 1).dot(w2[i].reshape(1, -1))
+          for i in xrange(len(w2))]) / len(w2)
+Q1 = np.matrix(Q1)
+Q2 = np.matrix(Q2)
+
+if len(w1) > 1:
+    while True:
+        try:
+            L1 = np.linalg.cholesky(Q1)
+            break
+        except np.linalg.linalg.LinAlgError:
+            mev = min(np.linalg.eig(Q1)[0]).real
+            assert (mev < 0)
+            print('Q1', mev)
+            if -mev < 1e-6:
+                mev = -1e-6
+            Q1 += (-mev * 1.000001) * np.identity(Q1.shape[0])
+
+    while True:
+        try:
+            L2 = np.linalg.cholesky(Q2)
+            break
+        except np.linalg.linalg.LinAlgError:
+            mev = min(np.linalg.eig(Q2)[0]).real
+            assert (mev < 0)
+            print('Q2', mev)
+            if -mev < 1e-6:
+                mev = -1e-6
+            Q2 += (-mev * 1.000001) * np.identity(Q2.shape[0])
+else:
+    L1 = np.identity(Q1.shape[0])
+    L2 = np.identity(Q2.shape[0])
+
+L1 = np.linalg.inv(L1) * gamma
+L2 = np.linalg.inv(L2) * gamma
+
+assert (L1.shape == L2.shape)
+assert (L1.shape[0] == L1.shape[1])
+
+f = open(args.output, 'wb')
+
+f.write(struct.pack('I', L1.shape[0]))
+f.write(struct.pack('I', L1.shape[1]))
+
+for i in xrange(L1.shape[0]):
+    for j in xrange(L1.shape[1]):
+        f.write(struct.pack('f', L1[i, j]))
+
+for i in xrange(L2.shape[0]):
+    for j in xrange(L2.shape[1]):
+        f.write(struct.pack('f', L2[i, j]))
+
+b1 = L1.dot(w1mean.reshape(-1, 1))
+b2 = L2.dot(w2mean.reshape(-1, 1))
+
+assert (L1.shape[0] == b1.shape[0])
+
+for i in xrange(b1.shape[0]):
+    f.write(struct.pack('f', b1[i, 0]))
+
+for i in xrange(b2.shape[0]):
+    f.write(struct.pack('f', b2[i, 0]))
+
+f.close()

modules/optflow/src/pcaflow.cpp

@@ -0,0 +1,526 @@
+/*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 "opencv2/ximgproc/edge_filter.hpp"
+#include "precomp.hpp"
+
+/* Disable "from double to float" and "from size_t to int" warnings.
+ * Fixing these would make the code look ugly by introducing explicit cast all around.
+ * Here these warning are pointless anyway.
+ */
+#ifdef _MSC_VER
+#pragma warning( disable : 4305 4244 4267 4838 )
+#endif
+#ifdef __clang__
+#pragma clang diagnostic ignored "-Wshorten-64-to-32"
+#endif
+
+namespace cv
+{
+namespace optflow
+{
+namespace
+{
+
+#ifndef M_SQRT2
+const float M_SQRT2 = 1.41421356237309504880;
+#endif
+
+template <typename T> inline int mathSign( T val ) { return ( T( 0 ) < val ) - ( val < T( 0 ) ); }
+
+/* Stable symmetric Householder reflection that gives c and s such that
+ *   [ c  s ][a] = [d],
+ *   [ s -c ][b]   [0]
+ *
+ * Output:
+ *   c -- cosine(theta), where theta is the implicit angle of rotation
+ *        (counter-clockwise) in a plane-rotation
+ *   s -- sine(theta)
+ *   r -- two-norm of [a; b]
+ */
+inline void symOrtho( double a, double b, double &c, double &s, double &r )
+{
+  if ( b == 0 )
+  {
+    c = mathSign( a );
+    s = 0;
+    r = std::abs( a );
+  }
+  else if ( a == 0 )
+  {
+    c = 0;
+    s = mathSign( b );
+    r = std::abs( b );
+  }
+  else if ( std::abs( b ) > std::abs( a ) )
+  {
+    const double tau = a / b;
+    s = mathSign( b ) / std::sqrt( 1 + tau * tau );
+    c = s * tau;
+    r = b / s;
+  }
+  else
+  {
+    const double tau = b / a;
+    c = mathSign( a ) / std::sqrt( 1 + tau * tau );
+    s = c * tau;
+    r = a / c;
+  }
+}
+
+/* Iterative LSQR algorithm for solving least squares problems.
+ *
+ * [1] Paige, C. C. and M. A. Saunders,
+ * LSQR: An Algorithm for Sparse Linear Equations And Sparse Least Squares
+ * ACM Trans. Math. Soft., Vol.8, 1982, pp. 43-71.
+ *
+ * Solves the following problem:
+ *   argmin_x ||Ax - b|| + damp||x||
+ *
+ * Output:
+ *   x -- approximate solution
+ */
+void solveLSQR( const Mat &A, const Mat &b, OutputArray xOut, const double damp = 0.0, const unsigned iter_lim = 10 )
+{
+  const int n = A.size().width;
+  CV_Assert( A.size().height == b.size().height );
+  CV_Assert( A.type() == CV_32F );
+  CV_Assert( b.type() == CV_32F );
+  xOut.create( n, 1, CV_32F );
+
+  Mat v( n, 1, CV_32F, 0.0f );
+  Mat u = b;
+  Mat x = xOut.getMat();
+  x = Mat::zeros( x.size(), x.type() );
+  double alfa = 0;
+  double beta = cv::norm( u, NORM_L2 );
+  Mat w( n, 1, CV_32F, 0.0f );
+  const Mat AT = A.t();
+
+  if ( beta > 0 )
+  {
+    u *= 1 / beta;
+    v = AT * u;
+    alfa = cv::norm( v, NORM_L2 );
+  }
+
+  if ( alfa > 0 )
+  {
+    v *= 1 / alfa;
+    w = v.clone();
+  }
+
+  double rhobar = alfa;
+  double phibar = beta;
+  if ( alfa * beta == 0 )
+    return;
+
+  for ( unsigned itn = 0; itn < iter_lim; ++itn )
+  {
+    u *= -alfa;
+    u += A * v;
+    beta = cv::norm( u, NORM_L2 );
+
+    if ( beta > 0 )
+    {
+      u *= 1 / beta;
+      v *= -beta;
+      v += AT * u;
+      alfa = cv::norm( v, NORM_L2 );
+      if ( alfa > 0 )
+        v *= 1 / alfa;
+    }
+
+    double rhobar1 = sqrt( rhobar * rhobar + damp * damp );
+    double cs1 = rhobar / rhobar1;
+    phibar = cs1 * phibar;
+
+    double cs, sn, rho;
+    symOrtho( rhobar1, beta, cs, sn, rho );
+
+    double theta = sn * alfa;
+    rhobar = -cs * alfa;
+    double phi = cs * phibar;
+    phibar = sn * phibar;
+
+    double t1 = phi / rho;
+    double t2 = -theta / rho;
+
+    x += t1 * w;
+    w *= t2;
+    w += v;
+  }
+}
+
+inline void _cpu_fillDCTSampledPoints( float *row, const Point2f &p, const Size &basisSize, const Size &size )
+{
+  for ( int n1 = 0; n1 < basisSize.width; ++n1 )
+    for ( int n2 = 0; n2 < basisSize.height; ++n2 )
+      row[n1 * basisSize.height + n2] =
+        cosf( ( n1 * CV_PI / size.width ) * ( p.x + 0.5 ) ) * cosf( ( n2 * CV_PI / size.height ) * ( p.y + 0.5 ) );
+}
+
+ocl::ProgramSource _ocl_fillDCTSampledPointsSource(
+  "__kernel void fillDCTSampledPoints(__global const uchar* features, int fstep, int foff, __global "
+  "uchar* A, int Astep, int Aoff, int fs, int bsw, int bsh, int sw, int sh) {"
+  "const int i = get_global_id(0);"
+  "const int n1 = get_global_id(1);"
+  "const int n2 = get_global_id(2);"
+  "if (i >= fs || n1 >= bsw || n2 >= bsh) return;"
+  "__global const float2* f = (__global const float2*)(features + (fstep * i + foff));"
+  "__global float* a = (__global float*)(A + (Astep * i + Aoff + (n1 * bsh + n2) * sizeof(float)));"
+  "const float2 p = f[0];"
+  "const float pi = 3.14159265358979323846;"
+  "a[0] = cos((n1 * pi / sw) * (p.x + 0.5)) * cos((n2 * pi / sh) * (p.y + 0.5));"
+  "}" );
+
+void applyCLAHE( UMat &img, float claheClip )
+{
+  Ptr<CLAHE> clahe = createCLAHE();
+  clahe->setClipLimit( claheClip );
+  clahe->apply( img, img );
+}
+
+void reduceToFlow( const Mat &w1, const Mat &w2, Mat &flow, const Size &basisSize )
+{
+  const Size size = flow.size();
+  Mat flowX( size, CV_32F, 0.0f );
+  Mat flowY( size, CV_32F, 0.0f );
+
+  const float mult = sqrt( size.area() ) * 0.5;
+
+  for ( int i = 0; i < basisSize.width; ++i )
+    for ( int j = 0; j < basisSize.height; ++j )
+    {
+      flowX.at<float>( j, i ) = w1.at<float>( i * basisSize.height + j ) * mult;
+      flowY.at<float>( j, i ) = w2.at<float>( i * basisSize.height + j ) * mult;
+    }
+  for ( int i = 0; i < basisSize.height; ++i )
+  {
+    flowX.at<float>( i, 0 ) *= M_SQRT2;
+    flowY.at<float>( i, 0 ) *= M_SQRT2;
+  }
+  for ( int i = 0; i < basisSize.width; ++i )
+  {
+    flowX.at<float>( 0, i ) *= M_SQRT2;
+    flowY.at<float>( 0, i ) *= M_SQRT2;
+  }
+
+  dct( flowX, flowX, DCT_INVERSE );
+  dct( flowY, flowY, DCT_INVERSE );
+  for ( int i = 0; i < size.height; ++i )
+    for ( int j = 0; j < size.width; ++j )
+      flow.at<Point2f>( i, j ) = Point2f( flowX.at<float>( i, j ), flowY.at<float>( i, j ) );
+}
+}
+
+void OpticalFlowPCAFlow::findSparseFeatures( UMat &from, UMat &to, std::vector<Point2f> &features,
+                                             std::vector<Point2f> &predictedFeatures ) const
+{
+  Size size = from.size();
+  const unsigned maxFeatures = size.area() * sparseRate;
+  goodFeaturesToTrack( from, features, maxFeatures * retainedCornersFraction, 0.005, 3 );
+
+  // Add points along the grid if not enough features
+  if ( maxFeatures > features.size() )
+  {
+    const unsigned missingPoints = maxFeatures - features.size();
+    const unsigned blockSize = sqrt( (float)size.area() / missingPoints );
+    for ( int x = blockSize / 2; x < size.width; x += blockSize )
+      for ( int y = blockSize / 2; y < size.height; y += blockSize )
+        features.push_back( Point2f( x, y ) );
+  }
+  std::vector<uchar> predictedStatus;
+  std::vector<float> predictedError;
+  calcOpticalFlowPyrLK( from, to, features, predictedFeatures, predictedStatus, predictedError );
+
+  size_t j = 0;
+  for ( size_t i = 0; i < features.size(); ++i )
+  {
+    if ( predictedStatus[i] )
+    {
+      features[j] = features[i];
+      predictedFeatures[j] = predictedFeatures[i];
+      ++j;
+    }
+  }
+  features.resize( j );
+  predictedFeatures.resize( j );
+}
+
+void OpticalFlowPCAFlow::removeOcclusions( UMat &from, UMat &to, std::vector<Point2f> &features,
+                                           std::vector<Point2f> &predictedFeatures ) const
+{
+  std::vector<uchar> predictedStatus;
+  std::vector<float> predictedError;
+  std::vector<Point2f> backwardFeatures;
+  calcOpticalFlowPyrLK( to, from, predictedFeatures, backwardFeatures, predictedStatus, predictedError );
+
+  size_t j = 0;
+  const float threshold = occlusionsThreshold * sqrt( from.size().area() );
+  for ( size_t i = 0; i < predictedFeatures.size(); ++i )
+  {
+    if ( predictedStatus[i] )
+    {
+      Point2f flowDiff = features[i] - backwardFeatures[i];
+      if ( flowDiff.dot( flowDiff ) <= threshold )
+      {
+        features[j] = features[i];
+        predictedFeatures[j] = predictedFeatures[i];
+        ++j;
+      }
+    }
+  }
+  features.resize( j );
+  predictedFeatures.resize( j );
+}
+
+void OpticalFlowPCAFlow::getSystem( OutputArray AOut, OutputArray b1Out, OutputArray b2Out,
+                                    const std::vector<Point2f> &features, const std::vector<Point2f> &predictedFeatures,
+                                    const Size size )
+{
+  AOut.create( features.size(), basisSize.area(), CV_32F );
+  b1Out.create( features.size(), 1, CV_32F );
+  b2Out.create( features.size(), 1, CV_32F );
+  if ( useOpenCL )
+  {
+    UMat A = AOut.getUMat();
+    Mat b1 = b1Out.getMat();
+    Mat b2 = b2Out.getMat();
+
+    ocl::Kernel kernel( "fillDCTSampledPoints", _ocl_fillDCTSampledPointsSource );
+    size_t globSize[] = {features.size(), basisSize.width, basisSize.height};
+    kernel
+      .args( cv::ocl::KernelArg::ReadOnlyNoSize( Mat( features ).getUMat( ACCESS_READ ) ),
+             cv::ocl::KernelArg::WriteOnlyNoSize( A ), (int)features.size(), (int)basisSize.width,
+             (int)basisSize.height, (int)size.width, (int)size.height )
+      .run( 3, globSize, 0, true );
+
+    for ( size_t i = 0; i < features.size(); ++i )
+    {
+      const Point2f flow = predictedFeatures[i] - features[i];
+      b1.at<float>( i ) = flow.x;
+      b2.at<float>( i ) = flow.y;
+    }
+  }
+  else
+  {
+    Mat A = AOut.getMat();
+    Mat b1 = b1Out.getMat();
+    Mat b2 = b2Out.getMat();
+
+    for ( size_t i = 0; i < features.size(); ++i )
+    {
+      _cpu_fillDCTSampledPoints( A.ptr<float>( i ), features[i], basisSize, size );
+      const Point2f flow = predictedFeatures[i] - features[i];
+      b1.at<float>( i ) = flow.x;
+      b2.at<float>( i ) = flow.y;
+    }
+  }
+}
+
+void OpticalFlowPCAFlow::getSystem( OutputArray A1Out, OutputArray A2Out, OutputArray b1Out, OutputArray b2Out,
+                                    const std::vector<Point2f> &features, const std::vector<Point2f> &predictedFeatures,
+                                    const Size size )
+{
+  CV_Assert( prior->getBasisSize() == basisSize.area() );
+
+  A1Out.create( features.size() + prior->getPadding(), basisSize.area(), CV_32F );
+  A2Out.create( features.size() + prior->getPadding(), basisSize.area(), CV_32F );
+  b1Out.create( features.size() + prior->getPadding(), 1, CV_32F );
+  b2Out.create( features.size() + prior->getPadding(), 1, CV_32F );
+
+  if ( useOpenCL )
+  {
+    UMat A = A1Out.getUMat();
+    Mat b1 = b1Out.getMat();
+    Mat b2 = b2Out.getMat();
+
+    ocl::Kernel kernel( "fillDCTSampledPoints", _ocl_fillDCTSampledPointsSource );
+    size_t globSize[] = {features.size(), basisSize.width, basisSize.height};
+    kernel
+      .args( cv::ocl::KernelArg::ReadOnlyNoSize( Mat( features ).getUMat( ACCESS_READ ) ),
+             cv::ocl::KernelArg::WriteOnlyNoSize( A ), (int)features.size(), (int)basisSize.width,
+             (int)basisSize.height, (int)size.width, (int)size.height )
+      .run( 3, globSize, 0, true );
+
+    for ( size_t i = 0; i < features.size(); ++i )
+    {
+      const Point2f flow = predictedFeatures[i] - features[i];
+      b1.at<float>( i ) = flow.x;
+      b2.at<float>( i ) = flow.y;
+    }
+  }
+  else
+  {
+    Mat A1 = A1Out.getMat();
+    Mat b1 = b1Out.getMat();
+    Mat b2 = b2Out.getMat();
+
+    for ( size_t i = 0; i < features.size(); ++i )
+    {
+      _cpu_fillDCTSampledPoints( A1.ptr<float>( i ), features[i], basisSize, size );
+      const Point2f flow = predictedFeatures[i] - features[i];
+      b1.at<float>( i ) = flow.x;
+      b2.at<float>( i ) = flow.y;
+    }
+  }
+
+  Mat A1 = A1Out.getMat();
+  Mat A2 = A2Out.getMat();
+  Mat b1 = b1Out.getMat();
+  Mat b2 = b2Out.getMat();
+
+  memcpy( A2.ptr<float>(), A1.ptr<float>(), features.size() * basisSize.area() * sizeof( float ) );
+  prior->fillConstraints( A1.ptr<float>( features.size(), 0 ), A2.ptr<float>( features.size(), 0 ),
+                          b1.ptr<float>( features.size(), 0 ), b2.ptr<float>( features.size(), 0 ) );
+}
+
+void OpticalFlowPCAFlow::calc( InputArray I0, InputArray I1, InputOutputArray flowOut )
+{
+  const Size size = I0.size();
+  CV_Assert( size == I1.size() );
+
+  UMat from, to;
+  if ( I0.channels() == 3 )
+  {
+    cvtColor( I0, from, COLOR_BGR2GRAY );
+    from.convertTo( from, CV_8U );
+  }
+  else
+  {
+    I0.getMat().convertTo( from, CV_8U );
+  }
+  if ( I1.channels() == 3 )
+  {
+    cvtColor( I1, to, COLOR_BGR2GRAY );
+    to.convertTo( to, CV_8U );
+  }
+  else
+  {
+    I1.getMat().convertTo( to, CV_8U );
+  }
+
+  CV_Assert( from.channels() == 1 );
+  CV_Assert( to.channels() == 1 );
+
+  const Mat fromOrig = from.getMat( ACCESS_READ ).clone();
+  useOpenCL = flowOut.isUMat() && ocl::useOpenCL();
+
+  applyCLAHE( from, claheClip );
+  applyCLAHE( to, claheClip );
+
+  std::vector<Point2f> features, predictedFeatures;
+  findSparseFeatures( from, to, features, predictedFeatures );
+  removeOcclusions( from, to, features, predictedFeatures );
+
+  flowOut.create( size, CV_32FC2 );
+  Mat flow = flowOut.getMat();
+
+  Mat w1, w2;
+  if ( prior.get() )
+  {
+    Mat A1, A2, b1, b2;
+    getSystem( A1, A2, b1, b2, features, predictedFeatures, size );
+    solveLSQR( A1, b1, w1, dampingFactor * size.area() );
+    solveLSQR( A2, b2, w2, dampingFactor * size.area() );
+  }
+  else
+  {
+    Mat A, b1, b2;
+    getSystem( A, b1, b2, features, predictedFeatures, size );
+    solveLSQR( A, b1, w1, dampingFactor * size.area() );
+    solveLSQR( A, b2, w2, dampingFactor * size.area() );
+  }
+  Mat flowSmall( ( size / 8 ) * 2, CV_32FC2 );
+  reduceToFlow( w1, w2, flowSmall, basisSize );
+  resize( flowSmall, flow, size, 0, 0, INTER_LINEAR );
+  ximgproc::fastGlobalSmootherFilter( fromOrig, flow, flow, 500, 2 );
+}
+
+OpticalFlowPCAFlow::OpticalFlowPCAFlow( Ptr<const PCAPrior> _prior, const Size _basisSize, float _sparseRate,
+                                        float _retainedCornersFraction, float _occlusionsThreshold,
+                                        float _dampingFactor, float _claheClip )
+    : prior( _prior ), basisSize( _basisSize ), sparseRate( _sparseRate ),
+      retainedCornersFraction( _retainedCornersFraction ), occlusionsThreshold( _occlusionsThreshold ),
+      dampingFactor( _dampingFactor ), claheClip( _claheClip ), useOpenCL( false )
+{
+  CV_Assert( sparseRate > 0 && sparseRate <= 0.1 );
+  CV_Assert( retainedCornersFraction >= 0 && retainedCornersFraction <= 1.0 );
+  CV_Assert( occlusionsThreshold > 0 );
+}
+
+void OpticalFlowPCAFlow::collectGarbage() {}
+
+Ptr<DenseOpticalFlow> createOptFlow_PCAFlow() { return makePtr<OpticalFlowPCAFlow>(); }
+
+PCAPrior::PCAPrior( const char *pathToPrior )
+{
+  FILE *f = fopen( pathToPrior, "rb" );
+  CV_Assert( f );
+
+  unsigned n = 0, m = 0;
+  CV_Assert( fread( &n, sizeof( n ), 1, f ) == 1 );
+  CV_Assert( fread( &m, sizeof( m ), 1, f ) == 1 );
+
+  L1.create( n, m, CV_32F );
+  L2.create( n, m, CV_32F );
+  c1.create( n, 1, CV_32F );
+  c2.create( n, 1, CV_32F );
+
+  CV_Assert( fread( L1.ptr<float>(), n * m * sizeof( float ), 1, f ) == 1 );
+  CV_Assert( fread( L2.ptr<float>(), n * m * sizeof( float ), 1, f ) == 1 );
+  CV_Assert( fread( c1.ptr<float>(), n * sizeof( float ), 1, f ) == 1 );
+  CV_Assert( fread( c2.ptr<float>(), n * sizeof( float ), 1, f ) == 1 );
+
+  fclose( f );
+}
+
+void PCAPrior::fillConstraints( float *A1, float *A2, float *b1, float *b2 ) const
+{
+  memcpy( A1, L1.ptr<float>(), L1.size().area() * sizeof( float ) );
+  memcpy( A2, L2.ptr<float>(), L2.size().area() * sizeof( float ) );
+  memcpy( b1, c1.ptr<float>(), c1.size().area() * sizeof( float ) );
+  memcpy( b2, c2.ptr<float>(), c2.size().area() * sizeof( float ) );
+}
+}
+}

modules/optflow/src/sparse_matching_gpc.cpp

@@ -0,0 +1,333 @@
+/*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 "opencv2/core/core_c.h"
+#include "opencv2/highgui.hpp"
+#include "precomp.hpp"
+
+namespace cv
+{
+namespace optflow
+{
+namespace
+{
+
+const int patchRadius = 10;
+const double thresholdMagnitudeFrac = 0.6666666666;
+const int globalIters = 3;
+const int localIters = 500;
+const int minNumberOfSamples = 2;
+//const bool debugOutput = true;
+
+struct Magnitude
+{
+  float val;
+  int i;
+  int j;
+
+  Magnitude( float _val, int _i, int _j ) : val( _val ), i( _i ), j( _j ) {}
+  Magnitude() {}
+
+  bool operator<( const Magnitude &m ) const { return val > m.val; }
+};
+
+struct PartitionPredicate1
+{
+  Vec< double, GPCPatchDescriptor::nFeatures > coef;
+  double rhs;
+
+  PartitionPredicate1( const Vec< double, GPCPatchDescriptor::nFeatures > &_coef, double _rhs ) : coef( _coef ), rhs( _rhs ) {}
+
+  bool operator()( const GPCPatchSample &sample ) const
+  {
+    const bool direction1 = ( coef.dot( sample.first.feature ) < rhs );
+    const bool direction2 = ( coef.dot( sample.second.feature ) < rhs );
+    return direction1 == false && direction1 == direction2;
+  }
+};
+
+struct PartitionPredicate2
+{
+  Vec< double, GPCPatchDescriptor::nFeatures > coef;
+  double rhs;
+
+  PartitionPredicate2( const Vec< double, GPCPatchDescriptor::nFeatures > &_coef, double _rhs ) : coef( _coef ), rhs( _rhs ) {}
+
+  bool operator()( const GPCPatchSample &sample ) const
+  {
+    const bool direction1 = ( coef.dot( sample.first.feature ) < rhs );
+    const bool direction2 = ( coef.dot( sample.second.feature ) < rhs );
+    return direction1 != direction2;
+  }
+};
+
+float normL2Sqr( const Vec2f &v ) { return v[0] * v[0] + v[1] * v[1]; }
+
+bool checkBounds( int i, int j, Size sz )
+{
+  return i >= patchRadius && j >= patchRadius && i + patchRadius < sz.height && j + patchRadius < sz.width;
+}
+
+void getTrainingSamples( const Mat &from, const Mat &to, const Mat &gt, GPCSamplesVector &samples )
+{
+  const Size sz = gt.size();
+  std::vector< Magnitude > mag;
+
+  for ( int i = patchRadius; i + patchRadius < sz.height; ++i )
+    for ( int j = patchRadius; j + patchRadius < sz.width; ++j )
+      mag.push_back( Magnitude( normL2Sqr( gt.at< Vec2f >( i, j ) ), i, j ) );
+
+  size_t n = size_t(mag.size() * thresholdMagnitudeFrac); // As suggested in the paper, we discard part of the training samples
+                                                          // with a small displacement and train to better distinguish hard pairs.
+  std::nth_element( mag.begin(), mag.begin() + n, mag.end() );
+  mag.resize( n );
+  std::random_shuffle( mag.begin(), mag.end() );
+  n /= patchRadius;
+  mag.resize( n );
+
+  Mat fromCh[3], toCh[3];
+  split( from, fromCh );
+  split( to, toCh );
+
+  for ( size_t k = 0; k < n; ++k )
+  {
+    int i0 = mag[k].i;
+    int j0 = mag[k].j;
+    int i1 = i0 + cvRound( gt.at< Vec2f >( i0, j0 )[1] );
+    int j1 = j0 + cvRound( gt.at< Vec2f >( i0, j0 )[0] );
+    if ( checkBounds( i1, j1, sz ) )
+      samples.push_back( std::make_pair( GPCPatchDescriptor( fromCh, i0, j0 ), GPCPatchDescriptor( toCh, i1, j1 ) ) );
+  }
+}
+
+/* Sample random number from Cauchy distribution. */
+double getRandomCauchyScalar()
+{
+  static RNG rng;
+  return tan( rng.uniform( -1.54, 1.54 ) ); // I intentionally used the value slightly less than PI/2 to enforce strictly
+                                            // zero probability for large numbers. Resulting PDF for Cauchy has
+                                            // truncated "tails".
+}
+
+/* Sample random vector from Cauchy distribution (pointwise, i.e. vector whose components are independent random
+ * variables from Cauchy distribution) */
+void getRandomCauchyVector( Vec< double, GPCPatchDescriptor::nFeatures > &v )
+{
+  for ( unsigned i = 0; i < GPCPatchDescriptor::nFeatures; ++i )
+    v[i] = getRandomCauchyScalar();
+}
+}
+
+GPCPatchDescriptor::GPCPatchDescriptor( const Mat *imgCh, int i, int j )
+{
+  Rect roi( j - patchRadius, i - patchRadius, 2 * patchRadius, 2 * patchRadius );
+  Mat freqDomain;
+  dct( imgCh[0]( roi ), freqDomain );
+
+  feature[0] = freqDomain.at< float >( 0, 0 );
+  feature[1] = freqDomain.at< float >( 0, 1 );
+  feature[2] = freqDomain.at< float >( 0, 2 );
+  feature[3] = freqDomain.at< float >( 0, 3 );
+
+  feature[4] = freqDomain.at< float >( 1, 0 );
+  feature[5] = freqDomain.at< float >( 1, 1 );
+  feature[6] = freqDomain.at< float >( 1, 2 );
+  feature[7] = freqDomain.at< float >( 1, 3 );
+
+  feature[8] = freqDomain.at< float >( 2, 0 );
+  feature[9] = freqDomain.at< float >( 2, 1 );
+  feature[10] = freqDomain.at< float >( 2, 2 );
+  feature[11] = freqDomain.at< float >( 2, 3 );
+
+  feature[12] = freqDomain.at< float >( 3, 0 );
+  feature[13] = freqDomain.at< float >( 3, 1 );
+  feature[14] = freqDomain.at< float >( 3, 2 );
+  feature[15] = freqDomain.at< float >( 3, 3 );
+
+  feature[16] = cv::sum( imgCh[1]( roi ) )[0] / ( 2 * patchRadius );
+  feature[17] = cv::sum( imgCh[2]( roi ) )[0] / ( 2 * patchRadius );
+}
+
+bool GPCTree::trainNode( size_t nodeId, SIter begin, SIter end, unsigned depth )
+{
+  if ( std::distance( begin, end ) < minNumberOfSamples )
+    return false;
+
+  if ( nodeId >= nodes.size() )
+    nodes.resize( nodeId + 1 );
+
+  Node &node = nodes[nodeId];
+
+  // Select the best hyperplane
+  unsigned globalBestScore = 0;
+  std::vector< double > values;
+
+  for ( int j = 0; j < globalIters; ++j )
+  { // Global search step
+    Vec< double, GPCPatchDescriptor::nFeatures > coef;
+    unsigned localBestScore = 0;
+    getRandomCauchyVector( coef );
+
+    for ( int i = 0; i < localIters; ++i )
+    { // Local search step
+      double randomModification = getRandomCauchyScalar();
+      const int pos = i % GPCPatchDescriptor::nFeatures;
+      std::swap( coef[pos], randomModification );
+      values.clear();
+
+      for ( SIter iter = begin; iter != end; ++iter )
+      {
+        values.push_back( coef.dot( iter->first.feature ) );
+        values.push_back( coef.dot( iter->second.feature ) );
+      }
+
+      std::nth_element( values.begin(), values.begin() + values.size() / 2, values.end() );
+      const double median = values[values.size() / 2];
+      unsigned correct = 0;
+
+      for ( SIter iter = begin; iter != end; ++iter )
+      {
+        const bool direction = ( coef.dot( iter->first.feature ) < median );
+        if ( direction == ( coef.dot( iter->second.feature ) < median ) )
+          ++correct;
+      }
+
+      if ( correct > localBestScore )
+        localBestScore = correct;
+      else
+        coef[pos] = randomModification;
+
+      if ( correct > globalBestScore )
+      {
+        globalBestScore = correct;
+        node.coef = coef;
+        node.rhs = median;
+
+        /*if ( debugOutput )
+        {
+          printf( "[%u] Updating weights: correct %.2f (%u/%ld)\n", depth, double( correct ) / std::distance( begin, end ), correct,
+                  std::distance( begin, end ) );
+          for ( unsigned k = 0; k < GPCPatchDescriptor::nFeatures; ++k )
+            printf( "%.3f ", coef[k] );
+          printf( "\n" );
+        }*/
+      }
+    }
+  }
+  // Partition vector with samples according to the hyperplane in QuickSort-like manner.
+  // Unlike QuickSort, we need to partition it into 3 parts (left subtree samples; undefined samples; right subtree
+  // samples), so we call it two times.
+  SIter leftEnd = std::partition( begin, end, PartitionPredicate1( node.coef, node.rhs ) ); // Separate left subtree samples from others.
+  SIter rightBegin =
+    std::partition( leftEnd, end, PartitionPredicate2( node.coef, node.rhs ) ); // Separate undefined samples from right subtree samples.
+
+  node.left = ( trainNode( nodeId * 2 + 1, begin, leftEnd, depth + 1 ) ) ? unsigned(nodeId * 2 + 1) : 0;
+  node.right = ( trainNode( nodeId * 2 + 2, rightBegin, end, depth + 1 ) ) ? unsigned(nodeId * 2 + 2) : 0;
+
+  return true;
+}
+
+void GPCTree::train( GPCSamplesVector &samples )
+{
+  nodes.reserve( samples.size() * 2 - 1 ); // set upper bound for the possible number of nodes so all subsequent resize() will be no-op
+  trainNode( 0, samples.begin(), samples.end(), 0 );
+}
+
+void GPCTree::write( FileStorage &fs ) const
+{
+  if ( nodes.empty() )
+    CV_Error( CV_StsBadArg, "Tree have not been trained" );
+  fs << "nodes" << nodes;
+}
+
+void GPCTree::read( const FileNode &fn ) { fn["nodes"] >> nodes; }
+
+Ptr< GPCTrainingSamples > GPCTrainingSamples::create( const std::vector< String > &imagesFrom, const std::vector< String > &imagesTo,
+                                                      const std::vector< String > &gt )
+{
+  CV_Assert( imagesFrom.size() == imagesTo.size() );
+  CV_Assert( imagesFrom.size() == gt.size() );
+
+  Ptr< GPCTrainingSamples > ts = makePtr< GPCTrainingSamples >();
+  for ( size_t i = 0; i < imagesFrom.size(); ++i )
+  {
+    Mat from = imread( imagesFrom[i] );
+    Mat to = imread( imagesTo[i] );
+    Mat gtFlow = readOpticalFlow( gt[i] );
+
+    CV_Assert( from.size == to.size );
+    CV_Assert( from.size == gtFlow.size );
+    CV_Assert( from.channels() == 3 );
+    CV_Assert( to.channels() == 3 );
+
+    from.convertTo( from, CV_32FC3 );
+    to.convertTo( to, CV_32FC3 );
+    cvtColor( from, from, COLOR_BGR2YCrCb );
+    cvtColor( to, to, COLOR_BGR2YCrCb );
+
+    getTrainingSamples( from, to, gtFlow, ts->samples );
+  }
+
+  return ts;
+}
+
+} // namespace optflow
+
+void write( FileStorage &fs, const String &name, const optflow::GPCTree::Node &node )
+{
+  cv::internal::WriteStructContext ws( fs, name, CV_NODE_SEQ + CV_NODE_FLOW );
+  for ( unsigned i = 0; i < optflow::GPCPatchDescriptor::nFeatures; ++i )
+    write( fs, node.coef[i] );
+  write( fs, node.rhs );
+  write( fs, (int)node.left );
+  write( fs, (int)node.right );
+}
+
+void read( const FileNode &fn, optflow::GPCTree::Node &node, optflow::GPCTree::Node )
+{
+  FileNodeIterator it = fn.begin();
+  for ( unsigned i = 0; i < optflow::GPCPatchDescriptor::nFeatures; ++i )
+    it >> node.coef[i];
+  it >> node.rhs >> (int &)node.left >> (int &)node.right;
+}
+
+} // namespace cv