diff --git a/modules/video/src/tvl1flow.cpp b/modules/video/src/tvl1flow.cpp index 426ff9f0a..aa2b9cd32 100644 --- a/modules/video/src/tvl1flow.cpp +++ b/modules/video/src/tvl1flow.cpp @@ -100,7 +100,7 @@ protected: double tau; double lambda; double theta; - double gamma; + double gamma; int nscales; int warps; double epsilon; @@ -111,7 +111,7 @@ protected: int medianFiltering; private: - void procOneScale(const Mat_& I0, const Mat_& I1, Mat_& u1, Mat_& u2); + void procOneScale(const Mat_& I0, const Mat_& I1, Mat_& u1, Mat_& u2, Mat_& u3); bool procOneScale_ocl(const UMat& I0, const UMat& I1, UMat& u1, UMat& u2); @@ -121,8 +121,8 @@ private: std::vector > I0s; std::vector > I1s; std::vector > u1s; - std::vector > u2s; - std::vector > u3s; + std::vector > u2s; + std::vector > u3s; Mat_ I1x_buf; Mat_ I1y_buf; @@ -138,25 +138,26 @@ private: Mat_ rho_c_buf; Mat_ v1_buf; - Mat_ v2_buf; + Mat_ v2_buf; + Mat_ v3_buf; Mat_ p11_buf; Mat_ p12_buf; Mat_ p21_buf; - Mat_ p22_buf; - Mat_ p31_buf; - Mat_ p32_buf; + Mat_ p22_buf; + Mat_ p31_buf; + Mat_ p32_buf; Mat_ div_p1_buf; - Mat_ div_p2_buf; - Mat_ div_p3_buf; + Mat_ div_p2_buf; + Mat_ div_p3_buf; Mat_ u1x_buf; Mat_ u1y_buf; Mat_ u2x_buf; - Mat_ u2y_buf; - Mat_ u3x_buf; - Mat_ u3y_buf; + Mat_ u2y_buf; + Mat_ u3x_buf; + Mat_ u3y_buf; } dm; struct dataUMat { @@ -350,7 +351,7 @@ OpticalFlowDual_TVL1::OpticalFlowDual_TVL1() nscales = 5; warps = 5; epsilon = 0.01; - gamma = 1.; + gamma = 0.01; innerIterations = 30; outerIterations = 10; useInitialFlow = false; @@ -375,15 +376,15 @@ void OpticalFlowDual_TVL1::calc(InputArray _I0, InputArray _I1, InputOutputArray dm.I0s.resize(nscales); dm.I1s.resize(nscales); dm.u1s.resize(nscales); - dm.u2s.resize(nscales); - dm.u3s.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()); - dm.u3s[0].create(I0.size()); + dm.u2s[0].create(I0.size()); + dm.u3s[0].create(I0.size()); if (useInitialFlow) { @@ -405,25 +406,26 @@ void OpticalFlowDual_TVL1::calc(InputArray _I0, InputArray _I1, InputOutputArray dm.rho_c_buf.create(I0.size()); dm.v1_buf.create(I0.size()); - dm.v2_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.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.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()); + 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) @@ -448,21 +450,21 @@ void OpticalFlowDual_TVL1::calc(InputArray _I0, InputArray _I1, InputOutputArray else { dm.u1s[s].create(dm.I0s[s].size()); - dm.u2s[s].create(dm.I0s[s].size()); + dm.u2s[s].create(dm.I0s[s].size()); } - dm.u3s[s].create(dm.I0s[s].size()); + 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)); - } - dm.u3s[nscales - 1].setTo(Scalar::all(0)); + } + 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) @@ -472,7 +474,7 @@ 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()); + resize(dm.u2s[s], dm.u2s[s - 1], dm.I0s[s - 1].size()); resize(dm.u3s[s], dm.u3s[s - 1], dm.I0s[s - 1].size()); // scale the optical flow with the appropriate zoom factor (don't scale u3!) @@ -888,10 +890,10 @@ void CalcGradRhoBody::operator() (const Range& range) const // compute the constant part of the rho function rhoRow[x] = (I1wRow[x] - I1wxRow[x] * u1Row[x] - I1wyRow[x] * u2Row[x] - I0Row[x]); - //It = I1wRow[x] - I0Row[x] - //(u - u0)*i_X = I1wxRow[x] * u1Row[x] - //(v - v0)*i_Y = I1wyRow[x] * u2Row[x] - // gamma * w = gamma * u3 + //It = I1wRow[x] - I0Row[x] + //(u - u0)*i_X = I1wxRow[x] * u1Row[x] + //(v - v0)*i_Y = I1wyRow[x] * u2Row[x] + // gamma * w = gamma * u3 } } } @@ -931,15 +933,15 @@ struct EstimateVBody : ParallelLoopBody Mat_ I1wx; Mat_ I1wy; Mat_ u1; - Mat_ u2; - Mat_ u3; + Mat_ u2; + Mat_ u3; Mat_ grad; Mat_ rho_c; mutable Mat_ v1; - mutable Mat_ v2; - mutable Mat_ v3; + mutable Mat_ v2; + mutable Mat_ v3; float l_t; - float gamma; + float gamma; }; void EstimateVBody::operator() (const Range& range) const @@ -949,14 +951,14 @@ void EstimateVBody::operator() (const Range& range) const const float* I1wxRow = I1wx[y]; const float* I1wyRow = I1wy[y]; const float* u1Row = u1[y]; - const float* u2Row = u2[y]; - const float* u3Row = u3[y]; + const float* u2Row = u2[y]; + const float* u3Row = u3[y]; const float* gradRow = grad[y]; const float* rhoRow = rho_c[y]; float* v1Row = v1[y]; - float* v2Row = v2[y]; - float* v3Row = v3[y]; + float* v2Row = v2[y]; + float* v3Row = v3[y]; for (int x = 0; x < I1wx.cols; ++x) { @@ -964,37 +966,37 @@ void EstimateVBody::operator() (const Range& range) const float d1 = 0.0f; float d2 = 0.0f; - float d3 = 0.0f; + float d3 = 0.0f; // add d3 for 3 cases if (rho < -l_t * gradRow[x]) { d1 = l_t * I1wxRow[x]; d2 = l_t * I1wyRow[x]; - d3 = l_t * gamma; + d3 = l_t * gamma; } else if (rho > l_t * gradRow[x]) { d1 = -l_t * I1wxRow[x]; - d2 = -l_t * I1wyRow[x]; - d3 = -l_t * gamma; + d2 = -l_t * I1wyRow[x]; + d3 = -l_t * gamma; } else if (gradRow[x] > std::numeric_limits::epsilon()) { float fi = -rho / gradRow[x]; d1 = fi * I1wxRow[x]; d2 = fi * I1wyRow[x]; - d3 = fi * gamma; + d3 = fi * gamma; } v1Row[x] = u1Row[x] + d1; - v2Row[x] = u2Row[x] + d2; - v3Row[x] = u3Row[x] + d3; + v2Row[x] = u2Row[x] + d2; + v3Row[x] = u3Row[x] + d3; } } } void estimateV(const Mat_& I1wx, const Mat_& I1wy, const Mat_& u1, const Mat_& u2, const Mat_& u3, const Mat_& grad, const Mat_& rho_c, - Mat_& v1, Mat_& v2, Mat_& v3, float l_t, float gamma) + Mat_& v1, Mat_& v2, Mat_& v3, float l_t, float gamma) { CV_DbgAssert( I1wy.size() == I1wx.size() ); CV_DbgAssert( u1.size() == I1wx.size() ); @@ -1011,15 +1013,15 @@ void estimateV(const Mat_& I1wx, const Mat_& I1wy, const Mat_& I1wx, const Mat_& I1wy, const Mat_& v1, const Mat_& v2, const Mat_& div_p1, const Mat_& div_p2, Mat_& u1, Mat_& u2, float theta) +float estimateU(const Mat_& v1, const Mat_& v2, const Mat_& v3, + const Mat_& div_p1, const Mat_& div_p2, const Mat_& div_p3, + Mat_& u1, Mat_& u2, Mat_& u3, + float theta, float gamma) { CV_DbgAssert( v2.size() == v1.size() ); + CV_DbgAssert( v3.size() == v1.size() ); CV_DbgAssert( div_p1.size() == v1.size() ); CV_DbgAssert( div_p2.size() == v1.size() ); + CV_DbgAssert( div_p3.size() == v1.size() ); CV_DbgAssert( u1.size() == v1.size() ); CV_DbgAssert( u2.size() == v1.size() ); + CV_DbgAssert( u3.size() == v1.size() ); float error = 0.0f; for (int y = 0; y < v1.rows; ++y) { const float* v1Row = v1[y]; - const float* v2Row = v2[y]; + const float* v2Row = v2[y]; + const float* v3Row = v3[y]; const float* divP1Row = div_p1[y]; - const float* divP2Row = div_p2[y]; + const float* divP2Row = div_p2[y]; + const float* divP3Row = div_p3[y]; float* u1Row = u1[y]; - float* u2Row = u2[y]; + float* u2Row = u2[y]; + float* u3Row = u3[y]; for (int x = 0; x < v1.cols; ++x) { const float u1k = u1Row[x]; - const float u2k = u2Row[x]; + const float u2k = u2Row[x]; + const float u3k = u3Row[x]; u1Row[x] = v1Row[x] + theta * divP1Row[x]; - u2Row[x] = v2Row[x] + theta * divP2Row[x]; - - //u3 - - error += (u1Row[x] - u1k) * (u1Row[x] - u1k) + (u2Row[x] - u2k) * (u2Row[x] - u2k); + u2Row[x] = v2Row[x] + theta * divP2Row[x]; + u3Row[x] = v3Row[x] + theta * divP3Row[x]; + + error += (u1Row[x] - u1k) * (u1Row[x] - u1k) + (u2Row[x] - u2k) * (u2Row[x] - u2k) + (u3Row[x] - u3k) * (u3Row[x] - u3k); } } @@ -1073,11 +1084,15 @@ struct EstimateDualVariablesBody : ParallelLoopBody Mat_ u1x; Mat_ u1y; Mat_ u2x; - Mat_ u2y; + Mat_ u2y; + Mat_ u3x; + Mat_ u3y; mutable Mat_ p11; mutable Mat_ p12; mutable Mat_ p21; - mutable Mat_ p22; + mutable Mat_ p22; + mutable Mat_ p31; + mutable Mat_ p32; float taut; }; @@ -1088,50 +1103,71 @@ void EstimateDualVariablesBody::operator() (const Range& range) const const float* u1xRow = u1x[y]; const float* u1yRow = u1y[y]; const float* u2xRow = u2x[y]; - const float* u2yRow = u2y[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* p22Row = p22[y]; + float* p31Row = p31[y]; + float* p32Row = p32[y]; for (int x = 0; x < u1x.cols; ++x) { const float g1 = static_cast(hypot(u1xRow[x], u1yRow[x])); - const float g2 = static_cast(hypot(u2xRow[x], u2yRow[x])); + const float g2 = static_cast(hypot(u2xRow[x], u2yRow[x])); + const float g3 = static_cast(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; + p22Row[x] = (p22Row[x] + taut * u2yRow[x]) / ng2; + p31Row[x] = (p31Row[x] + taut * u3xRow[x]) / ng3; + p32Row[x] = (p32Row[x] + taut * u3yRow[x]) / ng3; } } } -void estimateDualVariables(const Mat_& u1x, const Mat_& u1y, const Mat_& u2x, const Mat_& u2y, - Mat_& p11, Mat_& p12, Mat_& p21, Mat_& p22, float taut) +void estimateDualVariables(const Mat_& u1x, const Mat_& u1y, + const Mat_& u2x, const Mat_& u2y, + const Mat_& u3x, const Mat_& u3y, + Mat_& p11, Mat_& p12, + Mat_& p21, Mat_& p22, + Mat_& p31, Mat_& p32, + float taut) { 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; body.u1x = u1x; body.u1y = u1y; body.u2x = u2x; - body.u2y = u2y; + body.u2y = u2y; + body.u3x = u3x; + body.u3y = u3y; body.p11 = p11; body.p12 = p12; body.p21 = p21; - body.p22 = p22; + body.p22 = p22; + body.p31 = p31; + body.p32 = p32; body.taut = taut; parallel_for_(Range(0, u1x.rows), body); @@ -1225,7 +1261,7 @@ bool OpticalFlowDual_TVL1::procOneScale_ocl(const UMat& I0, const UMat& I1, UMat return true; } -void OpticalFlowDual_TVL1::procOneScale(const Mat_& I0, const Mat_& I1, Mat_& u1, Mat_& u2) +void OpticalFlowDual_TVL1::procOneScale(const Mat_& I0, const Mat_& I1, Mat_& u1, Mat_& u2, Mat_& u3) { const float scaledEpsilon = static_cast(epsilon * epsilon * I0.size().area()); @@ -1249,31 +1285,32 @@ void OpticalFlowDual_TVL1::procOneScale(const Mat_& I0, const Mat_ Mat_ rho_c = dm.rho_c_buf(Rect(0, 0, I0.cols, I0.rows)); Mat_ v1 = dm.v1_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ v2 = dm.v2_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ v2 = dm.v2_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ v3 = dm.v3_buf(Rect(0, 0, I0.cols, I0.rows)); Mat_ p11 = dm.p11_buf(Rect(0, 0, I0.cols, I0.rows)); Mat_ p12 = dm.p12_buf(Rect(0, 0, I0.cols, I0.rows)); Mat_ p21 = dm.p21_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ p22 = dm.p22_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ p31 = dm.p31_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ p32 = dm.p32_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ p22 = dm.p22_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ p31 = dm.p31_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ 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)); - p31.setTo(Scalar::all(0)); - p32.setTo(Scalar::all(0)); + p22.setTo(Scalar::all(0)); + p31.setTo(Scalar::all(0)); + p32.setTo(Scalar::all(0)); Mat_ div_p1 = dm.div_p1_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ div_p2 = dm.div_p2_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ div_p3 = dm.div_p2_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ div_p2 = dm.div_p2_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ div_p3 = dm.div_p3_buf(Rect(0, 0, I0.cols, I0.rows)); Mat_ u1x = dm.u1x_buf(Rect(0, 0, I0.cols, I0.rows)); Mat_ u1y = dm.u1y_buf(Rect(0, 0, I0.cols, I0.rows)); Mat_ u2x = dm.u2x_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ u2y = dm.u2y_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ u3x = dm.u3x_buf(Rect(0, 0, I0.cols, I0.rows)); - Mat_ u3y = dm.u3y_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ u2y = dm.u2y_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ u3x = dm.u3x_buf(Rect(0, 0, I0.cols, I0.rows)); + Mat_ u3y = dm.u3y_buf(Rect(0, 0, I0.cols, I0.rows)); const float l_t = static_cast(lambda * theta); const float taut = static_cast(tau / theta); @@ -1285,7 +1322,7 @@ void OpticalFlowDual_TVL1::procOneScale(const Mat_& I0, const Mat_ 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 + //calculate I1(x+u0) and its gradient calcGradRho(I0, I1w, I1wx, I1wy, u1, u2, grad, rho_c); float error = std::numeric_limits::max(); @@ -1298,21 +1335,23 @@ void OpticalFlowDual_TVL1::procOneScale(const Mat_& I0, const Mat_ 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, 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); + divergence(p21, p22, div_p2); + 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(theta)); + error = estimateU(v1, v2, v3, div_p1, div_p2, div_p3, u1, u2, u3, static_cast(theta), gamma); // compute the gradient of the optical flow (Du1, Du2) forwardGradient(u1, u1x, u1y); - forwardGradient(u2, u2x, u2y); + forwardGradient(u2, u2x, u2y); + 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); } } } @@ -1402,6 +1441,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 Ali term"); obj.info()->addParam(obj, "useInitialFlow", obj.useInitialFlow)) } // namespace