diff --git a/modules/cudaoptflow/include/opencv2/cudaoptflow.hpp b/modules/cudaoptflow/include/opencv2/cudaoptflow.hpp
index 220bb3457..d07a834ef 100644
--- a/modules/cudaoptflow/include/opencv2/cudaoptflow.hpp
+++ b/modules/cudaoptflow/include/opencv2/cudaoptflow.hpp
@@ -210,6 +210,14 @@ public:
      * In theory, it should have a small value in order to maintain both parts in correspondence.
      * The method is stable for a large range of values of this parameter.
      */
+
+    double gamma;
+    /**
+    * parameter used for motion estimation. It adds a variable allowing for illumination variations
+    * Set this parameter to 1. if you have varying illumination.
+    * See: Chambolle et al, A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging
+    * Journal of Mathematical imaging and vision, may 2011 Vol 40 issue 1, pp 120-145
+    */
     double theta;
 
     /**
@@ -241,12 +249,13 @@ public:
     bool useInitialFlow;
 
 private:
-    void procOneScale(const GpuMat& I0, const GpuMat& I1, GpuMat& u1, GpuMat& u2);
+    void procOneScale(const GpuMat& I0, const GpuMat& I1, GpuMat& u1, GpuMat& u2, GpuMat& u3);
 
     std::vector<GpuMat> I0s;
     std::vector<GpuMat> I1s;
     std::vector<GpuMat> u1s;
     std::vector<GpuMat> u2s;
+    std::vector<GpuMat> u3s;
 
     GpuMat I1x_buf;
     GpuMat I1y_buf;
@@ -262,6 +271,8 @@ private:
     GpuMat p12_buf;
     GpuMat p21_buf;
     GpuMat p22_buf;
+    GpuMat p31_buf;
+    GpuMat p32_buf;
 
     GpuMat diff_buf;
     GpuMat norm_buf;
diff --git a/modules/cudaoptflow/perf/perf_optflow.cpp b/modules/cudaoptflow/perf/perf_optflow.cpp
index 71ab89508..d22eb7e60 100644
--- a/modules/cudaoptflow/perf/perf_optflow.cpp
+++ b/modules/cudaoptflow/perf/perf_optflow.cpp
@@ -54,7 +54,7 @@ typedef pair<string, string> pair_string;
 DEF_PARAM_TEST_1(ImagePair, pair_string);
 
 PERF_TEST_P(ImagePair, InterpolateFrames,
-            Values<pair_string>(make_pair("gpu/opticalflow/frame0.png", "gpu/opticalflow/frame1.png")))
+    Values<pair_string>(make_pair("gpu/opticalflow/frame0.png", "gpu/opticalflow/frame1.png")))
 {
     cv::Mat frame0 = readImage(GetParam().first, cv::IMREAD_GRAYSCALE);
     ASSERT_FALSE(frame0.empty());
@@ -73,7 +73,7 @@ PERF_TEST_P(ImagePair, InterpolateFrames,
         cv::cuda::GpuMat d_bu, d_bv;
 
         cv::cuda::BroxOpticalFlow d_flow(0.197f /*alpha*/, 50.0f /*gamma*/, 0.8f /*scale_factor*/,
-                                        10 /*inner_iterations*/, 77 /*outer_iterations*/, 10 /*solver_iterations*/);
+            10 /*inner_iterations*/, 77 /*outer_iterations*/, 10 /*solver_iterations*/);
 
         d_flow(d_frame0, d_frame1, d_fu, d_fv);
         d_flow(d_frame1, d_frame0, d_bu, d_bv);
@@ -378,7 +378,6 @@ PERF_TEST_P(ImagePair, OpticalFlowDual_TVL1,
         alg->set("medianFiltering", 1);
         alg->set("innerIterations", 1);
         alg->set("outerIterations", 300);
-
         TEST_CYCLE() alg->calc(frame0, frame1, flow);
 
         CPU_SANITY_CHECK(flow);
@@ -389,7 +388,7 @@ PERF_TEST_P(ImagePair, OpticalFlowDual_TVL1,
 // OpticalFlowBM
 
 PERF_TEST_P(ImagePair, OpticalFlowBM,
-            Values<pair_string>(make_pair("gpu/opticalflow/frame0.png", "gpu/opticalflow/frame1.png")))
+    Values<pair_string>(make_pair("gpu/opticalflow/frame0.png", "gpu/opticalflow/frame1.png")))
 {
     declare.time(400);
 
@@ -421,7 +420,7 @@ PERF_TEST_P(ImagePair, OpticalFlowBM,
 }
 
 PERF_TEST_P(ImagePair, DISABLED_FastOpticalFlowBM,
-            Values<pair_string>(make_pair("gpu/opticalflow/frame0.png", "gpu/opticalflow/frame1.png")))
+    Values<pair_string>(make_pair("gpu/opticalflow/frame0.png", "gpu/opticalflow/frame1.png")))
 {
     declare.time(400);
 
diff --git a/modules/cudaoptflow/src/cuda/tvl1flow.cu b/modules/cudaoptflow/src/cuda/tvl1flow.cu
index b85dee701..2b66c972b 100644
--- a/modules/cudaoptflow/src/cuda/tvl1flow.cu
+++ b/modules/cudaoptflow/src/cuda/tvl1flow.cu
@@ -209,9 +209,11 @@ namespace tvl1flow
 
     __global__ void estimateUKernel(const PtrStepSzf I1wx, const PtrStepf I1wy,
                               const PtrStepf grad, const PtrStepf rho_c,
-                              const PtrStepf p11, const PtrStepf p12, const PtrStepf p21, const PtrStepf p22,
-                              PtrStepf u1, PtrStepf u2, PtrStepf error,
-                              const float l_t, const float theta, const bool calcError)
+                              const PtrStepf p11, const PtrStepf p12,
+                              const PtrStepf p21, const PtrStepf p22,
+                              const PtrStepf p31, const PtrStepf p32,
+                              PtrStepf u1, PtrStepf u2, PtrStepf u3, PtrStepf error,
+                              const float l_t, const float theta, const float gamma, const bool calcError)
     {
         const int x = blockIdx.x * blockDim.x + threadIdx.x;
         const int y = blockIdx.y * blockDim.y + threadIdx.y;
@@ -224,46 +226,59 @@ namespace tvl1flow
         const float gradVal = grad(y, x);
         const float u1OldVal = u1(y, x);
         const float u2OldVal = u2(y, x);
+        const float u3OldVal = gamma ? u3(y, x) : 0;
 
-        const float rho = rho_c(y, x) + (I1wxVal * u1OldVal + I1wyVal * u2OldVal);
+        const float rho = rho_c(y, x) + (I1wxVal * u1OldVal + I1wyVal * u2OldVal + gamma * u3OldVal);
 
         // estimate the values of the variable (v1, v2) (thresholding operator TH)
 
         float d1 = 0.0f;
         float d2 = 0.0f;
+        float d3 = 0.0f;
 
         if (rho < -l_t * gradVal)
         {
             d1 = l_t * I1wxVal;
             d2 = l_t * I1wyVal;
+            if (gamma)
+                d3 = l_t * gamma;
         }
         else if (rho > l_t * gradVal)
         {
             d1 = -l_t * I1wxVal;
             d2 = -l_t * I1wyVal;
+            if (gamma)
+                d3 = -l_t * gamma;
         }
         else if (gradVal > numeric_limits<float>::epsilon())
         {
             const float fi = -rho / gradVal;
             d1 = fi * I1wxVal;
             d2 = fi * I1wyVal;
+            if (gamma)
+                d3 = fi * gamma;
         }
 
         const float v1 = u1OldVal + d1;
         const float v2 = u2OldVal + d2;
+        const float v3 = u3OldVal + d3;
 
         // compute the divergence of the dual variable (p1, p2)
 
         const float div_p1 = divergence(p11, p12, y, x);
         const float div_p2 = divergence(p21, p22, y, x);
+        const float div_p3 = gamma ? divergence(p31, p32, y, x) : 0;
 
         // estimate the values of the optical flow (u1, u2)
 
         const float u1NewVal = v1 + theta * div_p1;
         const float u2NewVal = v2 + theta * div_p2;
+        const float u3NewVal = gamma ? v3 + theta * div_p3 : 0;
 
         u1(y, x) = u1NewVal;
         u2(y, x) = u2NewVal;
+        if (gamma)
+            u3(y, x) = u3NewVal;
 
         if (calcError)
         {
@@ -275,14 +290,14 @@ namespace tvl1flow
 
     void estimateU(PtrStepSzf I1wx, PtrStepSzf I1wy,
                    PtrStepSzf grad, PtrStepSzf rho_c,
-                   PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22,
-                   PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf error,
-                   float l_t, float theta, bool calcError)
+                   PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32,
+                   PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf error,
+                   float l_t, float theta, float gamma, bool calcError)
     {
         const dim3 block(32, 8);
         const dim3 grid(divUp(I1wx.cols, block.x), divUp(I1wx.rows, block.y));
 
-        estimateUKernel<<<grid, block>>>(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, u1, u2, error, l_t, theta, calcError);
+        estimateUKernel<<<grid, block>>>(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, p31, p32, u1, u2, u3, error, l_t, theta, gamma, calcError);
         cudaSafeCall( cudaGetLastError() );
 
         cudaSafeCall( cudaDeviceSynchronize() );
@@ -294,7 +309,8 @@ namespace tvl1flow
 
 namespace tvl1flow
 {
-    __global__ void estimateDualVariablesKernel(const PtrStepSzf u1, const PtrStepf u2, PtrStepf p11, PtrStepf p12, PtrStepf p21, PtrStepf p22, const float taut)
+    __global__ void estimateDualVariablesKernel(const PtrStepSzf u1, const PtrStepf u2, const PtrStepSzf u3,
+                                                PtrStepf p11, PtrStepf p12, PtrStepf p21, PtrStepf p22, PtrStepf p31, PtrStepf p32, const float taut, const float gamma)
     {
         const int x = blockIdx.x * blockDim.x + threadIdx.x;
         const int y = blockIdx.y * blockDim.y + threadIdx.y;
@@ -308,24 +324,34 @@ namespace tvl1flow
         const float u2x = u2(y, ::min(x + 1, u1.cols - 1)) - u2(y, x);
         const float u2y = u2(::min(y + 1, u1.rows - 1), x) - u2(y, x);
 
+        const float u3x = gamma ? u3(y, ::min(x + 1, u1.cols - 1)) - u3(y, x) : 0;
+        const float u3y = gamma ? u3(::min(y + 1, u1.rows - 1), x) - u3(y, x) : 0;
+
         const float g1 = ::hypotf(u1x, u1y);
         const float g2 = ::hypotf(u2x, u2y);
+        const float g3 = gamma ? ::hypotf(u3x, u3y) : 0;
 
         const float ng1 = 1.0f + taut * g1;
         const float ng2 = 1.0f + taut * g2;
+        const float ng3 = gamma ? 1.0f + taut * g3 : 0;
 
         p11(y, x) = (p11(y, x) + taut * u1x) / ng1;
         p12(y, x) = (p12(y, x) + taut * u1y) / ng1;
         p21(y, x) = (p21(y, x) + taut * u2x) / ng2;
         p22(y, x) = (p22(y, x) + taut * u2y) / ng2;
+        if (gamma)
+        {
+            p31(y, x) = (p31(y, x) + taut * u3x) / ng3;
+            p32(y, x) = (p32(y, x) + taut * u3y) / ng3;
+        }
     }
 
-    void estimateDualVariables(PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, float taut)
+    void estimateDualVariables(PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32, float taut, float gamma)
     {
         const dim3 block(32, 8);
         const dim3 grid(divUp(u1.cols, block.x), divUp(u1.rows, block.y));
 
-        estimateDualVariablesKernel<<<grid, block>>>(u1, u2, p11, p12, p21, p22, taut);
+        estimateDualVariablesKernel<<<grid, block>>>(u1, u2, u3, p11, p12, p21, p22, p31, p32, taut, gamma);
         cudaSafeCall( cudaGetLastError() );
 
         cudaSafeCall( cudaDeviceSynchronize() );
diff --git a/modules/cudaoptflow/src/tvl1flow.cpp b/modules/cudaoptflow/src/tvl1flow.cpp
index 7b6882d9f..a54e212ba 100644
--- a/modules/cudaoptflow/src/tvl1flow.cpp
+++ b/modules/cudaoptflow/src/tvl1flow.cpp
@@ -64,6 +64,7 @@ cv::cuda::OpticalFlowDual_TVL1_CUDA::OpticalFlowDual_TVL1_CUDA()
     epsilon        = 0.01;
     iterations     = 300;
     scaleStep      = 0.8;
+    gamma           = 0.0;
     useInitialFlow = false;
 }
 
@@ -80,6 +81,7 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::operator ()(const GpuMat& I0, const Gp
     I1s.resize(nscales);
     u1s.resize(nscales);
     u2s.resize(nscales);
+    u3s.resize(nscales);
 
     I0.convertTo(I0s[0], CV_32F, I0.depth() == CV_8U ? 1.0 : 255.0);
     I1.convertTo(I1s[0], CV_32F, I1.depth() == CV_8U ? 1.0 : 255.0);
@@ -92,6 +94,8 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::operator ()(const GpuMat& I0, const Gp
 
     u1s[0] = flowx;
     u2s[0] = flowy;
+    if (gamma)
+        u3s[0].create(I0.size(), CV_32FC1);
 
     I1x_buf.create(I0.size(), CV_32FC1);
     I1y_buf.create(I0.size(), CV_32FC1);
@@ -107,7 +111,11 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::operator ()(const GpuMat& I0, const Gp
     p12_buf.create(I0.size(), CV_32FC1);
     p21_buf.create(I0.size(), CV_32FC1);
     p22_buf.create(I0.size(), CV_32FC1);
-
+    if (gamma)
+    {
+        p31_buf.create(I0.size(), CV_32FC1);
+        p32_buf.create(I0.size(), CV_32FC1);
+    }
     diff_buf.create(I0.size(), CV_32FC1);
 
     // create the scales
@@ -135,6 +143,8 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::operator ()(const GpuMat& I0, const Gp
             u1s[s].create(I0s[s].size(), CV_32FC1);
             u2s[s].create(I0s[s].size(), CV_32FC1);
         }
+        if (gamma)
+            u3s[s].create(I0s[s].size(), CV_32FC1);
     }
 
     if (!useInitialFlow)
@@ -142,12 +152,14 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::operator ()(const GpuMat& I0, const Gp
         u1s[nscales-1].setTo(Scalar::all(0));
         u2s[nscales-1].setTo(Scalar::all(0));
     }
+    if (gamma)
+        u3s[nscales - 1].setTo(Scalar::all(0));
 
     // pyramidal structure for computing the optical flow
     for (int s = nscales - 1; s >= 0; --s)
     {
         // compute the optical flow at the current scale
-        procOneScale(I0s[s], I1s[s], u1s[s], u2s[s]);
+        procOneScale(I0s[s], I1s[s], u1s[s], u2s[s], u3s[s]);
 
         // if this was the last scale, finish now
         if (s == 0)
@@ -158,6 +170,8 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::operator ()(const GpuMat& I0, const Gp
         // zoom the optical flow for the next finer scale
         cuda::resize(u1s[s], u1s[s - 1], I0s[s - 1].size());
         cuda::resize(u2s[s], u2s[s - 1], I0s[s - 1].size());
+        if (gamma)
+            cuda::resize(u3s[s], u3s[s - 1], I0s[s - 1].size());
 
         // scale the optical flow with the appropriate zoom factor
         cuda::multiply(u1s[s - 1], Scalar::all(1/scaleStep), u1s[s - 1]);
@@ -171,13 +185,13 @@ namespace tvl1flow
     void warpBackward(PtrStepSzf I0, PtrStepSzf I1, PtrStepSzf I1x, PtrStepSzf I1y, PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf I1w, PtrStepSzf I1wx, PtrStepSzf I1wy, PtrStepSzf grad, PtrStepSzf rho);
     void estimateU(PtrStepSzf I1wx, PtrStepSzf I1wy,
                    PtrStepSzf grad, PtrStepSzf rho_c,
-                   PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22,
-                   PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf error,
-                   float l_t, float theta, bool calcError);
-    void estimateDualVariables(PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, float taut);
+                   PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32,
+                   PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf error,
+                   float l_t, float theta, float gamma, bool calcError);
+    void estimateDualVariables(PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32, float taut, const float gamma);
 }
 
-void cv::cuda::OpticalFlowDual_TVL1_CUDA::procOneScale(const GpuMat& I0, const GpuMat& I1, GpuMat& u1, GpuMat& u2)
+void cv::cuda::OpticalFlowDual_TVL1_CUDA::procOneScale(const GpuMat& I0, const GpuMat& I1, GpuMat& u1, GpuMat& u2, GpuMat& u3)
 {
     using namespace tvl1flow;
 
@@ -203,10 +217,21 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::procOneScale(const GpuMat& I0, const G
     GpuMat p12 = p12_buf(Rect(0, 0, I0.cols, I0.rows));
     GpuMat p21 = p21_buf(Rect(0, 0, I0.cols, I0.rows));
     GpuMat p22 = p22_buf(Rect(0, 0, I0.cols, I0.rows));
+    GpuMat p31, p32;
+    if (gamma)
+    {
+        p31 = p31_buf(Rect(0, 0, I0.cols, I0.rows));
+        p32 = p32_buf(Rect(0, 0, I0.cols, I0.rows));
+    }
     p11.setTo(Scalar::all(0));
     p12.setTo(Scalar::all(0));
     p21.setTo(Scalar::all(0));
     p22.setTo(Scalar::all(0));
+    if (gamma)
+    {
+        p31.setTo(Scalar::all(0));
+        p32.setTo(Scalar::all(0));
+    }
 
     GpuMat diff = diff_buf(Rect(0, 0, I0.cols, I0.rows));
 
@@ -223,9 +248,8 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::procOneScale(const GpuMat& I0, const G
         {
             // some tweaks to make sum operation less frequently
             bool calcError = (epsilon > 0) && (n & 0x1) && (prevError < scaledEpsilon);
-
-            estimateU(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, u1, u2, diff, l_t, static_cast<float>(theta), calcError);
-
+            cv::Mat m1(u3);
+            estimateU(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, p31, p32, u1, u2, u3, diff, l_t, static_cast<float>(theta), gamma, calcError);
             if (calcError)
             {
                 error = cuda::sum(diff, norm_buf)[0];
@@ -237,7 +261,7 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::procOneScale(const GpuMat& I0, const G
                 prevError -= scaledEpsilon;
             }
 
-            estimateDualVariables(u1, u2, p11, p12, p21, p22, taut);
+            estimateDualVariables(u1, u2, u3, p11, p12, p21, p22, p31, p32, taut, gamma);
         }
     }
 }
@@ -248,6 +272,7 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::collectGarbage()
     I1s.clear();
     u1s.clear();
     u2s.clear();
+    u3s.clear();
 
     I1x_buf.release();
     I1y_buf.release();
@@ -263,7 +288,11 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::collectGarbage()
     p12_buf.release();
     p21_buf.release();
     p22_buf.release();
-
+    if (gamma)
+    {
+        p31_buf.release();
+        p32_buf.release();
+    }
     diff_buf.release();
     norm_buf.release();
 }
diff --git a/modules/cudaoptflow/test/test_optflow.cpp b/modules/cudaoptflow/test/test_optflow.cpp
index 1de40510d..2b976563b 100644
--- a/modules/cudaoptflow/test/test_optflow.cpp
+++ b/modules/cudaoptflow/test/test_optflow.cpp
@@ -360,6 +360,18 @@ CUDA_TEST_P(OpticalFlowDual_TVL1, Accuracy)
     alg->calc(frame0, frame1, flow);
     cv::Mat gold[2];
     cv::split(flow, gold);
+    cv::Mat mx(d_flowx);
+    cv::Mat my(d_flowx);
+
+    EXPECT_MAT_SIMILAR(gold[0], d_flowx, 4e-3);
+    EXPECT_MAT_SIMILAR(gold[1], d_flowy, 4e-3);
+    d_alg.gamma = 1;
+    alg->set("gamma", 1);
+    d_alg(loadMat(frame0, useRoi), loadMat(frame1, useRoi), d_flowx, d_flowy);
+    alg->calc(frame0, frame1, flow);
+    cv::split(flow, gold);
+    mx = cv::Mat(d_flowx);
+    my = cv::Mat(d_flowx);
 
     EXPECT_MAT_SIMILAR(gold[0], d_flowx, 4e-3);
     EXPECT_MAT_SIMILAR(gold[1], d_flowy, 4e-3);
diff --git a/modules/cudastereo/src/cuda/stereocsbp.cu b/modules/cudastereo/src/cuda/stereocsbp.cu
index dd535e8b2..394693240 100644
--- a/modules/cudastereo/src/cuda/stereocsbp.cu
+++ b/modules/cudastereo/src/cuda/stereocsbp.cu
@@ -61,7 +61,7 @@ namespace cv { namespace cuda { namespace device
         template <int channels> static float __device__ pixeldiff(const uchar* left, const uchar* right, float max_data_term);
         template<> __device__ __forceinline__ static float pixeldiff<1>(const uchar* left, const uchar* right, float max_data_term)
         {
-            return fmin( ::abs((int)*left - *right), max_data_term);
+            return fminf( ::abs((int)*left - *right), max_data_term);
         }
         template<> __device__ __forceinline__ static float pixeldiff<3>(const uchar* left, const uchar* right, float max_data_term)
         {
@@ -69,7 +69,7 @@ namespace cv { namespace cuda { namespace device
             float tg = 0.587f * ::abs((int)left[1] - right[1]);
             float tr = 0.299f * ::abs((int)left[2] - right[2]);
 
-            return fmin(tr + tg + tb, max_data_term);
+            return fminf(tr + tg + tb, max_data_term);
         }
         template<> __device__ __forceinline__ static float pixeldiff<4>(const uchar* left, const uchar* right, float max_data_term)
         {
@@ -80,7 +80,7 @@ namespace cv { namespace cuda { namespace device
             float tg = 0.587f * ::abs((int)l.y - r.y);
             float tr = 0.299f * ::abs((int)l.z - r.z);
 
-            return fmin(tr + tg + tb, max_data_term);
+            return fminf(tr + tg + tb, max_data_term);
         }
 
         template <typename T>
diff --git a/modules/video/src/tvl1flow.cpp b/modules/video/src/tvl1flow.cpp
index fec000dc4..8a865a1d2 100644
--- a/modules/video/src/tvl1flow.cpp
+++ b/modules/video/src/tvl1flow.cpp
@@ -100,6 +100,7 @@ protected:
     double tau;
     double lambda;
     double theta;
+    double gamma;
     int nscales;
     int warps;
     double epsilon;
@@ -110,7 +111,7 @@ protected:
     int medianFiltering;
 
 private:
-    void procOneScale(const Mat_<float>& I0, const Mat_<float>& I1, Mat_<float>& u1, Mat_<float>& u2);
+   void procOneScale(const Mat_<float>& I0, const Mat_<float>& I1, Mat_<float>& u1, Mat_<float>& u2, Mat_<float>& u3);
 
     bool procOneScale_ocl(const UMat& I0, const UMat& I1, UMat& u1, UMat& u2);
 
@@ -121,6 +122,7 @@ private:
         std::vector<Mat_<float> > I1s;
         std::vector<Mat_<float> > u1s;
         std::vector<Mat_<float> > u2s;
+        std::vector<Mat_<float> > u3s;
 
         Mat_<float> I1x_buf;
         Mat_<float> I1y_buf;
@@ -137,19 +139,25 @@ private:
 
         Mat_<float> v1_buf;
         Mat_<float> v2_buf;
+        Mat_<float> v3_buf;
 
         Mat_<float> p11_buf;
         Mat_<float> p12_buf;
         Mat_<float> p21_buf;
         Mat_<float> p22_buf;
+        Mat_<float> p31_buf;
+        Mat_<float> p32_buf;
 
         Mat_<float> div_p1_buf;
         Mat_<float> div_p2_buf;
+        Mat_<float> div_p3_buf;
 
         Mat_<float> u1x_buf;
         Mat_<float> u1y_buf;
         Mat_<float> u2x_buf;
         Mat_<float> u2y_buf;
+        Mat_<float> u3x_buf;
+        Mat_<float> u3y_buf;
     } dm;
     struct dataUMat
     {
@@ -343,6 +351,7 @@ OpticalFlowDual_TVL1::OpticalFlowDual_TVL1()
     nscales        = 5;
     warps          = 5;
     epsilon        = 0.01;
+    gamma          = 0.;
     innerIterations = 30;
     outerIterations = 10;
     useInitialFlow = false;
@@ -364,18 +373,20 @@ void OpticalFlowDual_TVL1::calc(InputArray _I0, InputArray _I1, InputOutputArray
     CV_Assert( I0.type() == I1.type() );
     CV_Assert( !useInitialFlow || (_flow.size() == I0.size() && _flow.type() == CV_32FC2) );
     CV_Assert( nscales > 0 );
-
+    bool use_gamma = gamma != 0;
     // allocate memory for the pyramid structure
     dm.I0s.resize(nscales);
     dm.I1s.resize(nscales);
     dm.u1s.resize(nscales);
     dm.u2s.resize(nscales);
+    dm.u3s.resize(nscales);
 
     I0.convertTo(dm.I0s[0], dm.I0s[0].depth(), I0.depth() == CV_8U ? 1.0 : 255.0);
     I1.convertTo(dm.I1s[0], dm.I1s[0].depth(), I1.depth() == CV_8U ? 1.0 : 255.0);
 
     dm.u1s[0].create(I0.size());
     dm.u2s[0].create(I0.size());
+    if (use_gamma) dm.u3s[0].create(I0.size());
 
     if (useInitialFlow)
     {
@@ -398,19 +409,25 @@ void OpticalFlowDual_TVL1::calc(InputArray _I0, InputArray _I1, InputOutputArray
 
     dm.v1_buf.create(I0.size());
     dm.v2_buf.create(I0.size());
+    dm.v3_buf.create(I0.size());
 
     dm.p11_buf.create(I0.size());
     dm.p12_buf.create(I0.size());
     dm.p21_buf.create(I0.size());
     dm.p22_buf.create(I0.size());
+    dm.p31_buf.create(I0.size());
+    dm.p32_buf.create(I0.size());
 
     dm.div_p1_buf.create(I0.size());
     dm.div_p2_buf.create(I0.size());
+    dm.div_p3_buf.create(I0.size());
 
     dm.u1x_buf.create(I0.size());
     dm.u1y_buf.create(I0.size());
     dm.u2x_buf.create(I0.size());
     dm.u2y_buf.create(I0.size());
+    dm.u3x_buf.create(I0.size());
+    dm.u3y_buf.create(I0.size());
 
     // create the scales
     for (int s = 1; s < nscales; ++s)
@@ -437,19 +454,19 @@ void OpticalFlowDual_TVL1::calc(InputArray _I0, InputArray _I1, InputOutputArray
             dm.u1s[s].create(dm.I0s[s].size());
             dm.u2s[s].create(dm.I0s[s].size());
         }
+        if (use_gamma) dm.u3s[s].create(dm.I0s[s].size());
     }
-
     if (!useInitialFlow)
     {
         dm.u1s[nscales - 1].setTo(Scalar::all(0));
         dm.u2s[nscales - 1].setTo(Scalar::all(0));
     }
-
+    if (use_gamma) dm.u3s[nscales - 1].setTo(Scalar::all(0));
     // pyramidal structure for computing the optical flow
     for (int s = nscales - 1; s >= 0; --s)
     {
         // compute the optical flow at the current scale
-        procOneScale(dm.I0s[s], dm.I1s[s], dm.u1s[s], dm.u2s[s]);
+        procOneScale(dm.I0s[s], dm.I1s[s], dm.u1s[s], dm.u2s[s], dm.u3s[s]);
 
         // if this was the last scale, finish now
         if (s == 0)
@@ -460,8 +477,9 @@ void OpticalFlowDual_TVL1::calc(InputArray _I0, InputArray _I1, InputOutputArray
         // zoom the optical flow for the next finer scale
         resize(dm.u1s[s], dm.u1s[s - 1], dm.I0s[s - 1].size());
         resize(dm.u2s[s], dm.u2s[s - 1], dm.I0s[s - 1].size());
+        if (use_gamma) resize(dm.u3s[s], dm.u3s[s - 1], dm.I0s[s - 1].size());
 
-        // scale the optical flow with the appropriate zoom factor
+        // scale the optical flow with the appropriate zoom factor (don't scale u3!)
         multiply(dm.u1s[s - 1], Scalar::all(1 / scaleStep), dm.u1s[s - 1]);
         multiply(dm.u2s[s - 1], Scalar::all(1 / scaleStep), dm.u2s[s - 1]);
     }
@@ -914,59 +932,69 @@ struct EstimateVBody : ParallelLoopBody
     Mat_<float> I1wy;
     Mat_<float> u1;
     Mat_<float> u2;
+    Mat_<float> u3;
     Mat_<float> grad;
     Mat_<float> rho_c;
     mutable Mat_<float> v1;
     mutable Mat_<float> v2;
+    mutable Mat_<float> v3;
     float l_t;
+    float gamma;
 };
 
 void EstimateVBody::operator() (const Range& range) const
 {
+    bool use_gamma = gamma != 0;
     for (int y = range.start; y < range.end; ++y)
     {
         const float* I1wxRow = I1wx[y];
         const float* I1wyRow = I1wy[y];
         const float* u1Row = u1[y];
         const float* u2Row = u2[y];
+        const float* u3Row = use_gamma?u3[y]:NULL;
         const float* gradRow = grad[y];
         const float* rhoRow = rho_c[y];
 
         float* v1Row = v1[y];
         float* v2Row = v2[y];
+        float* v3Row = use_gamma ? v3[y]:NULL;
 
         for (int x = 0; x < I1wx.cols; ++x)
         {
-            const float rho = rhoRow[x] + (I1wxRow[x] * u1Row[x] + I1wyRow[x] * u2Row[x]);
-
+            const float rho = use_gamma ? rhoRow[x] + (I1wxRow[x] * u1Row[x] + I1wyRow[x] * u2Row[x]) + gamma * u3Row[x] :
+                                          rhoRow[x] + (I1wxRow[x] * u1Row[x] + I1wyRow[x] * u2Row[x]);
             float d1 = 0.0f;
             float d2 = 0.0f;
-
+            float d3 = 0.0f;
             if (rho < -l_t * gradRow[x])
             {
                 d1 = l_t * I1wxRow[x];
                 d2 = l_t * I1wyRow[x];
+                if (use_gamma) d3 = l_t * gamma;
             }
             else if (rho > l_t * gradRow[x])
             {
                 d1 = -l_t * I1wxRow[x];
                 d2 = -l_t * I1wyRow[x];
+                if (use_gamma) d3 = -l_t * gamma;
             }
             else if (gradRow[x] > std::numeric_limits<float>::epsilon())
             {
                 float fi = -rho / gradRow[x];
                 d1 = fi * I1wxRow[x];
                 d2 = fi * I1wyRow[x];
+                if (use_gamma) d3 = fi * gamma;
             }
 
             v1Row[x] = u1Row[x] + d1;
             v2Row[x] = u2Row[x] + d2;
+            if (use_gamma) v3Row[x] = u3Row[x] + d3;
         }
     }
 }
 
-void estimateV(const Mat_<float>& I1wx, const Mat_<float>& I1wy, const Mat_<float>& u1, const Mat_<float>& u2, const Mat_<float>& grad, const Mat_<float>& rho_c,
-               Mat_<float>& v1, Mat_<float>& v2, float l_t)
+void estimateV(const Mat_<float>& I1wx, const Mat_<float>& I1wy, const Mat_<float>& u1, const Mat_<float>& u2, const Mat_<float>& u3, const Mat_<float>& grad, const Mat_<float>& rho_c,
+   Mat_<float>& v1, Mat_<float>& v2, Mat_<float>& v3, float l_t, float gamma)
 {
     CV_DbgAssert( I1wy.size() == I1wx.size() );
     CV_DbgAssert( u1.size() == I1wx.size() );
@@ -977,24 +1005,29 @@ void estimateV(const Mat_<float>& I1wx, const Mat_<float>& I1wy, const Mat_<floa
     CV_DbgAssert( v2.size() == I1wx.size() );
 
     EstimateVBody body;
-
+    bool use_gamma = gamma != 0;
     body.I1wx = I1wx;
     body.I1wy = I1wy;
     body.u1 = u1;
     body.u2 = u2;
+    if (use_gamma) body.u3 = u3;
     body.grad = grad;
     body.rho_c = rho_c;
     body.v1 = v1;
     body.v2 = v2;
+    if (use_gamma) body.v3 = v3;
     body.l_t = l_t;
-
+    body.gamma = gamma;
     parallel_for_(Range(0, I1wx.rows), body);
 }
 
 ////////////////////////////////////////////////////////////
 // estimateU
 
-float estimateU(const Mat_<float>& v1, const Mat_<float>& v2, const Mat_<float>& div_p1, const Mat_<float>& div_p2, Mat_<float>& u1, Mat_<float>& u2, float theta)
+float estimateU(const Mat_<float>& v1, const Mat_<float>& v2, const Mat_<float>& v3,
+            const Mat_<float>& div_p1, const Mat_<float>& div_p2, const Mat_<float>& div_p3,
+            Mat_<float>& u1, Mat_<float>& u2, Mat_<float>& u3,
+            float theta, float gamma)
 {
     CV_DbgAssert( v2.size() == v1.size() );
     CV_DbgAssert( div_p1.size() == v1.size() );
@@ -1003,25 +1036,32 @@ float estimateU(const Mat_<float>& v1, const Mat_<float>& v2, const Mat_<float>&
     CV_DbgAssert( u2.size() == v1.size() );
 
     float error = 0.0f;
+    bool use_gamma = gamma != 0;
     for (int y = 0; y < v1.rows; ++y)
     {
         const float* v1Row = v1[y];
         const float* v2Row = v2[y];
+        const float* v3Row = use_gamma?v3[y]:NULL;
         const float* divP1Row = div_p1[y];
         const float* divP2Row = div_p2[y];
+        const float* divP3Row = use_gamma?div_p3[y]:NULL;
 
         float* u1Row = u1[y];
         float* u2Row = u2[y];
+        float* u3Row = use_gamma?u3[y]:NULL;
+
 
         for (int x = 0; x < v1.cols; ++x)
         {
             const float u1k = u1Row[x];
             const float u2k = u2Row[x];
+            const float u3k = use_gamma?u3Row[x]:0;
 
             u1Row[x] = v1Row[x] + theta * divP1Row[x];
             u2Row[x] = v2Row[x] + theta * divP2Row[x];
-
-            error += (u1Row[x] - u1k) * (u1Row[x] - u1k) + (u2Row[x] - u2k) * (u2Row[x] - u2k);
+            if (use_gamma) u3Row[x] = v3Row[x] + theta * divP3Row[x];
+            error += use_gamma?(u1Row[x] - u1k) * (u1Row[x] - u1k) + (u2Row[x] - u2k) * (u2Row[x] - u2k) + (u3Row[x] - u3k) * (u3Row[x] - u3k):
+                               (u1Row[x] - u1k) * (u1Row[x] - u1k) + (u2Row[x] - u2k) * (u2Row[x] - u2k);
         }
     }
 
@@ -1039,11 +1079,16 @@ struct EstimateDualVariablesBody : ParallelLoopBody
     Mat_<float> u1y;
     Mat_<float> u2x;
     Mat_<float> u2y;
+    Mat_<float> u3x;
+    Mat_<float> u3y;
     mutable Mat_<float> p11;
     mutable Mat_<float> p12;
     mutable Mat_<float> p21;
     mutable Mat_<float> p22;
+    mutable Mat_<float> p31;
+    mutable Mat_<float> p32;
     float taut;
+    bool use_gamma;
 };
 
 void EstimateDualVariablesBody::operator() (const Range& range) const
@@ -1054,38 +1099,55 @@ void EstimateDualVariablesBody::operator() (const Range& range) const
         const float* u1yRow = u1y[y];
         const float* u2xRow = u2x[y];
         const float* u2yRow = u2y[y];
+        const float* u3xRow = u3x[y];
+        const float* u3yRow = u3y[y];
 
         float* p11Row = p11[y];
         float* p12Row = p12[y];
         float* p21Row = p21[y];
         float* p22Row = p22[y];
+        float* p31Row = p31[y];
+        float* p32Row = p32[y];
 
         for (int x = 0; x < u1x.cols; ++x)
         {
             const float g1 = static_cast<float>(hypot(u1xRow[x], u1yRow[x]));
             const float g2 = static_cast<float>(hypot(u2xRow[x], u2yRow[x]));
+            const float g3 = static_cast<float>(hypot(u3xRow[x], u3yRow[x]));
 
             const float ng1  = 1.0f + taut * g1;
-            const float ng2  = 1.0f + taut * g2;
+            const float ng2 =  1.0f + taut * g2;
+            const float ng3 = 1.0f + taut * g3;
 
             p11Row[x] = (p11Row[x] + taut * u1xRow[x]) / ng1;
             p12Row[x] = (p12Row[x] + taut * u1yRow[x]) / ng1;
             p21Row[x] = (p21Row[x] + taut * u2xRow[x]) / ng2;
             p22Row[x] = (p22Row[x] + taut * u2yRow[x]) / ng2;
+            if (use_gamma) p31Row[x] = (p31Row[x] + taut * u3xRow[x]) / ng3;
+            if (use_gamma) p32Row[x] = (p32Row[x] + taut * u3yRow[x]) / ng3;
         }
     }
 }
 
-void estimateDualVariables(const Mat_<float>& u1x, const Mat_<float>& u1y, const Mat_<float>& u2x, const Mat_<float>& u2y,
-                           Mat_<float>& p11, Mat_<float>& p12, Mat_<float>& p21, Mat_<float>& p22, float taut)
+void estimateDualVariables(const Mat_<float>& u1x, const Mat_<float>& u1y,
+                     const Mat_<float>& u2x, const Mat_<float>& u2y,
+                     const Mat_<float>& u3x, const Mat_<float>& u3y,
+                           Mat_<float>& p11, Mat_<float>& p12,
+                     Mat_<float>& p21, Mat_<float>& p22,
+                     Mat_<float>& p31, Mat_<float>& p32,
+                     float taut, bool use_gamma)
 {
     CV_DbgAssert( u1y.size() == u1x.size() );
     CV_DbgAssert( u2x.size() == u1x.size() );
+    CV_DbgAssert( u3x.size() == u1x.size() );
     CV_DbgAssert( u2y.size() == u1x.size() );
+    CV_DbgAssert( u3y.size() == u1x.size() );
     CV_DbgAssert( p11.size() == u1x.size() );
     CV_DbgAssert( p12.size() == u1x.size() );
     CV_DbgAssert( p21.size() == u1x.size() );
     CV_DbgAssert( p22.size() == u1x.size() );
+    CV_DbgAssert( p31.size() == u1x.size() );
+    CV_DbgAssert( p32.size() == u1x.size() );
 
     EstimateDualVariablesBody body;
 
@@ -1093,11 +1155,16 @@ void estimateDualVariables(const Mat_<float>& u1x, const Mat_<float>& u1y, const
     body.u1y = u1y;
     body.u2x = u2x;
     body.u2y = u2y;
+    body.u3x = u3x;
+    body.u3y = u3y;
     body.p11 = p11;
     body.p12 = p12;
     body.p21 = p21;
     body.p22 = p22;
+    body.p31 = p31;
+    body.p32 = p32;
     body.taut = taut;
+    body.use_gamma = use_gamma;
 
     parallel_for_(Range(0, u1x.rows), body);
 }
@@ -1190,7 +1257,7 @@ bool OpticalFlowDual_TVL1::procOneScale_ocl(const UMat& I0, const UMat& I1, UMat
     return true;
 }
 
-void OpticalFlowDual_TVL1::procOneScale(const Mat_<float>& I0, const Mat_<float>& I1, Mat_<float>& u1, Mat_<float>& u2)
+void OpticalFlowDual_TVL1::procOneScale(const Mat_<float>& I0, const Mat_<float>& I1, Mat_<float>& u1, Mat_<float>& u2, Mat_<float>& u3)
 {
     const float scaledEpsilon = static_cast<float>(epsilon * epsilon * I0.size().area());
 
@@ -1215,23 +1282,32 @@ void OpticalFlowDual_TVL1::procOneScale(const Mat_<float>& I0, const Mat_<float>
 
     Mat_<float> v1 = dm.v1_buf(Rect(0, 0, I0.cols, I0.rows));
     Mat_<float> v2 = dm.v2_buf(Rect(0, 0, I0.cols, I0.rows));
+    Mat_<float> v3 = dm.v3_buf(Rect(0, 0, I0.cols, I0.rows));
 
     Mat_<float> p11 = dm.p11_buf(Rect(0, 0, I0.cols, I0.rows));
     Mat_<float> p12 = dm.p12_buf(Rect(0, 0, I0.cols, I0.rows));
     Mat_<float> p21 = dm.p21_buf(Rect(0, 0, I0.cols, I0.rows));
     Mat_<float> p22 = dm.p22_buf(Rect(0, 0, I0.cols, I0.rows));
+    Mat_<float> p31 = dm.p31_buf(Rect(0, 0, I0.cols, I0.rows));
+    Mat_<float> p32 = dm.p32_buf(Rect(0, 0, I0.cols, I0.rows));
     p11.setTo(Scalar::all(0));
     p12.setTo(Scalar::all(0));
     p21.setTo(Scalar::all(0));
     p22.setTo(Scalar::all(0));
+    bool use_gamma = gamma != 0.;
+    if (use_gamma) p31.setTo(Scalar::all(0));
+    if (use_gamma) p32.setTo(Scalar::all(0));
 
     Mat_<float> div_p1 = dm.div_p1_buf(Rect(0, 0, I0.cols, I0.rows));
     Mat_<float> div_p2 = dm.div_p2_buf(Rect(0, 0, I0.cols, I0.rows));
+    Mat_<float> div_p3 = dm.div_p3_buf(Rect(0, 0, I0.cols, I0.rows));
 
     Mat_<float> u1x = dm.u1x_buf(Rect(0, 0, I0.cols, I0.rows));
     Mat_<float> u1y = dm.u1y_buf(Rect(0, 0, I0.cols, I0.rows));
     Mat_<float> u2x = dm.u2x_buf(Rect(0, 0, I0.cols, I0.rows));
     Mat_<float> u2y = dm.u2y_buf(Rect(0, 0, I0.cols, I0.rows));
+    Mat_<float> u3x = dm.u3x_buf(Rect(0, 0, I0.cols, I0.rows));
+    Mat_<float> u3y = dm.u3y_buf(Rect(0, 0, I0.cols, I0.rows));
 
     const float l_t = static_cast<float>(lambda * theta);
     const float taut = static_cast<float>(tau / theta);
@@ -1243,7 +1319,7 @@ void OpticalFlowDual_TVL1::procOneScale(const Mat_<float>& I0, const Mat_<float>
         remap(I1, I1w, flowMap1, flowMap2, INTER_CUBIC);
         remap(I1x, I1wx, flowMap1, flowMap2, INTER_CUBIC);
         remap(I1y, I1wy, flowMap1, flowMap2, INTER_CUBIC);
-
+        //calculate I1(x+u0) and its gradient
         calcGradRho(I0, I1w, I1wx, I1wy, u1, u2, grad, rho_c);
 
         float error = std::numeric_limits<float>::max();
@@ -1256,21 +1332,23 @@ void OpticalFlowDual_TVL1::procOneScale(const Mat_<float>& I0, const Mat_<float>
             for (int n_inner = 0; error > scaledEpsilon && n_inner < innerIterations; ++n_inner)
             {
                 // estimate the values of the variable (v1, v2) (thresholding operator TH)
-                estimateV(I1wx, I1wy, u1, u2, grad, rho_c, v1, v2, l_t);
+                estimateV(I1wx, I1wy, u1, u2, u3, grad, rho_c, v1, v2, v3, l_t, static_cast<float>(gamma));
 
-                // compute the divergence of the dual variable (p1, p2)
+                // compute the divergence of the dual variable (p1, p2, p3)
                 divergence(p11, p12, div_p1);
                 divergence(p21, p22, div_p2);
+                if (use_gamma) divergence(p31, p32, div_p3);
 
                 // estimate the values of the optical flow (u1, u2)
-                error = estimateU(v1, v2, div_p1, div_p2, u1, u2, static_cast<float>(theta));
+                error = estimateU(v1, v2, v3, div_p1, div_p2, div_p3, u1, u2, u3, static_cast<float>(theta), static_cast<float>(gamma));
 
                 // compute the gradient of the optical flow (Du1, Du2)
                 forwardGradient(u1, u1x, u1y);
                 forwardGradient(u2, u2x, u2y);
+                if (use_gamma) forwardGradient(u3, u3x, u3y);
 
-                // estimate the values of the dual variable (p1, p2)
-                estimateDualVariables(u1x, u1y, u2x, u2y, p11, p12, p21, p22, taut);
+                // estimate the values of the dual variable (p1, p2, p3)
+                estimateDualVariables(u1x, u1y, u2x, u2y, u3x, u3y, p11, p12, p21, p22, p31, p32, taut, use_gamma);
             }
         }
     }
@@ -1360,6 +1438,8 @@ CV_INIT_ALGORITHM(OpticalFlowDual_TVL1, "DenseOpticalFlow.DualTVL1",
                                        "inner iterations (between outlier filtering) used in the numerical scheme");
                   obj.info()->addParam(obj, "outerIterations", obj.outerIterations, false, 0, 0,
                                        "outer iterations (number of inner loops) used in the numerical scheme");
+                  obj.info()->addParam(obj, "gamma", obj.gamma, false, 0, 0,
+                                       "coefficient for additional illumination variation term");
                   obj.info()->addParam(obj, "useInitialFlow", obj.useInitialFlow))
 
 } // namespace
diff --git a/modules/video/test/test_tvl1optflow.cpp b/modules/video/test/test_tvl1optflow.cpp
index 4772f0fb6..b829c7c89 100644
--- a/modules/video/test/test_tvl1optflow.cpp
+++ b/modules/video/test/test_tvl1optflow.cpp
@@ -166,7 +166,7 @@ TEST(Video_calcOpticalFlowDual_TVL1, Regression)
     ASSERT_EQ(gold.rows, flow.rows);
     ASSERT_EQ(gold.cols, flow.cols);
 
-    const double err = calcRMSE(gold, flow);
+    double err = calcRMSE(gold, flow);
     EXPECT_LE(err, MAX_RMSE);
 #endif
 }