From a3a09cf4d1eaf0bb46192886d5793b4acf2eedf7 Mon Sep 17 00:00:00 2001 From: Vladislav Vinogradov Date: Thu, 14 Feb 2013 13:26:47 +0400 Subject: [PATCH] refactored OpticalFlowDual_TVL1: * added DenseOpticalFlow interface * moved OpticalFlowDual_TVL1 to src folder --- modules/gpu/perf/perf_video.cpp | 6 +- modules/gpu/test/test_optflow.cpp | 4 +- .../doc/motion_analysis_and_object_tracking.rst | 18 +- modules/video/include/opencv2/video/tracking.hpp | 101 +-- modules/video/perf/perf_tvl1optflow.cpp | 7 +- modules/video/src/tvl1flow.cpp | 956 +++++++++++---------- modules/video/test/test_tvl1optflow.cpp | 4 +- samples/cpp/tvl1_optical_flow.cpp | 4 +- 8 files changed, 542 insertions(+), 558 deletions(-) diff --git a/modules/gpu/perf/perf_video.cpp b/modules/gpu/perf/perf_video.cpp index bf2fd99..83213a1 100644 --- a/modules/gpu/perf/perf_video.cpp +++ b/modules/gpu/perf/perf_video.cpp @@ -431,13 +431,13 @@ PERF_TEST_P(ImagePair, Video_OpticalFlowDual_TVL1, { cv::Mat flow; - cv::OpticalFlowDual_TVL1 alg; + cv::Ptr alg = cv::createOptFlow_DualTVL1(); - alg(frame0, frame1, flow); + alg->calc(frame0, frame1, flow); TEST_CYCLE() { - alg(frame0, frame1, flow); + alg->calc(frame0, frame1, flow); } CPU_SANITY_CHECK(flow); diff --git a/modules/gpu/test/test_optflow.cpp b/modules/gpu/test/test_optflow.cpp index c93ebbe..a97516f 100644 --- a/modules/gpu/test/test_optflow.cpp +++ b/modules/gpu/test/test_optflow.cpp @@ -431,9 +431,9 @@ GPU_TEST_P(OpticalFlowDual_TVL1, Accuracy) cv::gpu::GpuMat d_flowy = createMat(frame0.size(), CV_32FC1, useRoi); d_alg(loadMat(frame0, useRoi), loadMat(frame1, useRoi), d_flowx, d_flowy); - cv::OpticalFlowDual_TVL1 alg; + cv::Ptr alg = cv::createOptFlow_DualTVL1(); cv::Mat flow; - alg(frame0, frame1, flow); + alg->calc(frame0, frame1, flow); cv::Mat gold[2]; cv::split(flow, gold); diff --git a/modules/video/doc/motion_analysis_and_object_tracking.rst b/modules/video/doc/motion_analysis_and_object_tracking.rst index 2372373..2c5dcb1 100644 --- a/modules/video/doc/motion_analysis_and_object_tracking.rst +++ b/modules/video/doc/motion_analysis_and_object_tracking.rst @@ -643,11 +643,11 @@ See [Tao2012]_. And site of project - http://graphics.berkeley.edu/papers/Tao-SA -OpticalFlowDual_TVL1 --------------------- +createOptFlow_DualTVL1 +---------------------- "Dual TV L1" Optical Flow Algorithm. -.. ocv:class:: OpticalFlowDual_TVL12 +.. ocv:function:: Ptr createOptFlow_DualTVL1() The class implements the "Dual TV L1" optical flow algorithm described in [Zach2007]_ and [Javier2012]_ . @@ -685,11 +685,11 @@ Here are important members of the class that control the algorithm, which you ca -OpticalFlowDual_TVL1::operator() --------------------------------- +DenseOpticalFlow::calc +-------------------------- Calculates an optical flow. -.. ocv:function:: void OpticalFlowDual_TVL1::operator ()(InputArray I0, InputArray I1, InputOutputArray flow) +.. ocv:function:: void DenseOpticalFlow::calc(InputArray I0, InputArray I1, InputOutputArray flow) :param prev: first 8-bit single-channel input image. @@ -699,11 +699,11 @@ Calculates an optical flow. -OpticalFlowDual_TVL1::collectGarbage ------------------------------------- +DenseOpticalFlow::collectGarbage +-------------------------------- Releases all inner buffers. -.. ocv:function:: void OpticalFlowDual_TVL1::collectGarbage() +.. ocv:function:: void DenseOpticalFlow::collectGarbage() diff --git a/modules/video/include/opencv2/video/tracking.hpp b/modules/video/include/opencv2/video/tracking.hpp index deef20b..48e84c8 100644 --- a/modules/video/include/opencv2/video/tracking.hpp +++ b/modules/video/include/opencv2/video/tracking.hpp @@ -352,104 +352,19 @@ CV_EXPORTS_W void calcOpticalFlowSF(Mat& from, double upscale_sigma_color, double speed_up_thr); +class CV_EXPORTS DenseOpticalFlow : public Algorithm +{ +public: + virtual void calc(InputArray I0, InputArray I1, InputOutputArray flow) = 0; + virtual void collectGarbage() = 0; +}; + // Implementation of the Zach, Pock and Bischof Dual TV-L1 Optical Flow method // // see reference: // [1] C. Zach, T. Pock and H. Bischof, "A Duality Based Approach for Realtime TV-L1 Optical Flow". // [2] Javier Sanchez, Enric Meinhardt-Llopis and Gabriele Facciolo. "TV-L1 Optical Flow Estimation". -class CV_EXPORTS OpticalFlowDual_TVL1 -{ -public: - OpticalFlowDual_TVL1(); - - void operator ()(InputArray I0, InputArray I1, InputOutputArray flow); - - void collectGarbage(); - - /** - * Time step of the numerical scheme. - */ - double tau; - - /** - * Weight parameter for the data term, attachment parameter. - * This is the most relevant parameter, which determines the smoothness of the output. - * The smaller this parameter is, the smoother the solutions we obtain. - * It depends on the range of motions of the images, so its value should be adapted to each image sequence. - */ - double lambda; - - /** - * Weight parameter for (u - v)^2, tightness parameter. - * It serves as a link between the attachment and the regularization terms. - * 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 theta; - - /** - * Number of scales used to create the pyramid of images. - */ - int nscales; - - /** - * Number of warpings per scale. - * Represents the number of times that I1(x+u0) and grad( I1(x+u0) ) are computed per scale. - * This is a parameter that assures the stability of the method. - * It also affects the running time, so it is a compromise between speed and accuracy. - */ - int warps; - - /** - * Stopping criterion threshold used in the numerical scheme, which is a trade-off between precision and running time. - * A small value will yield more accurate solutions at the expense of a slower convergence. - */ - double epsilon; - - /** - * Stopping criterion iterations number used in the numerical scheme. - */ - int iterations; - - bool useInitialFlow; - -private: - void procOneScale(const Mat_& I0, const Mat_& I1, Mat_& u1, Mat_& u2); - - std::vector > I0s; - std::vector > I1s; - std::vector > u1s; - std::vector > u2s; - - Mat_ I1x_buf; - Mat_ I1y_buf; - - Mat_ flowMap1_buf; - Mat_ flowMap2_buf; - - Mat_ I1w_buf; - Mat_ I1wx_buf; - Mat_ I1wy_buf; - - Mat_ grad_buf; - Mat_ rho_c_buf; - - Mat_ v1_buf; - Mat_ v2_buf; - - Mat_ p11_buf; - Mat_ p12_buf; - Mat_ p21_buf; - Mat_ p22_buf; - - Mat_ div_p1_buf; - Mat_ div_p2_buf; - - Mat_ u1x_buf; - Mat_ u1y_buf; - Mat_ u2x_buf; - Mat_ u2y_buf; -}; +CV_EXPORTS Ptr createOptFlow_DualTVL1(); } diff --git a/modules/video/perf/perf_tvl1optflow.cpp b/modules/video/perf/perf_tvl1optflow.cpp index 2014130..ad90915 100644 --- a/modules/video/perf/perf_tvl1optflow.cpp +++ b/modules/video/perf/perf_tvl1optflow.cpp @@ -22,12 +22,9 @@ PERF_TEST_P(ImagePair, OpticalFlowDual_TVL1, testing::Values(impair("cv/optflow/ Mat flow; - OpticalFlowDual_TVL1 tvl1; + Ptr tvl1 = createOptFlow_DualTVL1(); - TEST_CYCLE() - { - tvl1(frame1, frame2, flow); - } + TEST_CYCLE_N(10) tvl1->calc(frame1, frame2, flow); SANITY_CHECK(flow, 0.5); } diff --git a/modules/video/src/tvl1flow.cpp b/modules/video/src/tvl1flow.cpp index 2f8034a..bff1d7e 100644 --- a/modules/video/src/tvl1flow.cpp +++ b/modules/video/src/tvl1flow.cpp @@ -77,7 +77,67 @@ using namespace std; using namespace cv; -cv::OpticalFlowDual_TVL1::OpticalFlowDual_TVL1() +namespace { + +class OpticalFlowDual_TVL1 : public DenseOpticalFlow +{ +public: + OpticalFlowDual_TVL1(); + + void calc(InputArray I0, InputArray I1, InputOutputArray flow); + void collectGarbage(); + + AlgorithmInfo* info() const; + +protected: + double tau; + double lambda; + double theta; + int nscales; + int warps; + double epsilon; + int iterations; + bool useInitialFlow; + +private: + void procOneScale(const Mat_& I0, const Mat_& I1, Mat_& u1, Mat_& u2); + + std::vector > I0s; + std::vector > I1s; + std::vector > u1s; + std::vector > u2s; + + Mat_ I1x_buf; + Mat_ I1y_buf; + + Mat_ flowMap1_buf; + Mat_ flowMap2_buf; + + Mat_ I1w_buf; + Mat_ I1wx_buf; + Mat_ I1wy_buf; + + Mat_ grad_buf; + Mat_ rho_c_buf; + + Mat_ v1_buf; + Mat_ v2_buf; + + Mat_ p11_buf; + Mat_ p12_buf; + Mat_ p21_buf; + Mat_ p22_buf; + + Mat_ div_p1_buf; + Mat_ div_p2_buf; + + Mat_ u1x_buf; + Mat_ u1y_buf; + Mat_ u2x_buf; + Mat_ u2y_buf; +}; + +OpticalFlowDual_TVL1::OpticalFlowDual_TVL1() { tau = 0.25; lambda = 0.15; @@ -89,7 +149,7 @@ cv::OpticalFlowDual_TVL1::OpticalFlowDual_TVL1() useInitialFlow = false; } -void cv::OpticalFlowDual_TVL1::operator ()(InputArray _I0, InputArray _I1, InputOutputArray _flow) +void OpticalFlowDual_TVL1::calc(InputArray _I0, InputArray _I1, InputOutputArray _flow) { Mat I0 = _I0.getMat(); Mat I1 = _I1.getMat(); @@ -195,542 +255,539 @@ void cv::OpticalFlowDual_TVL1::operator ()(InputArray _I0, InputArray _I1, Input merge(uxy, 2, _flow); } -namespace -{ - //////////////////////////////////////////////////////////// - // buildFlowMap +//////////////////////////////////////////////////////////// +// buildFlowMap - struct BuildFlowMapBody : ParallelLoopBody - { - void operator() (const Range& range) const; +struct BuildFlowMapBody : ParallelLoopBody +{ + void operator() (const Range& range) const; - Mat_ u1; - Mat_ u2; - mutable Mat_ map1; - mutable Mat_ map2; - }; + Mat_ u1; + Mat_ u2; + mutable Mat_ map1; + mutable Mat_ map2; +}; - void BuildFlowMapBody::operator() (const Range& range) const +void BuildFlowMapBody::operator() (const Range& range) const +{ + for (int y = range.start; y < range.end; ++y) { - for (int y = range.start; y < range.end; ++y) - { - const float* u1Row = u1[y]; - const float* u2Row = u2[y]; + const float* u1Row = u1[y]; + const float* u2Row = u2[y]; - float* map1Row = map1[y]; - float* map2Row = map2[y]; + float* map1Row = map1[y]; + float* map2Row = map2[y]; - for (int x = 0; x < u1.cols; ++x) - { - map1Row[x] = x + u1Row[x]; - map2Row[x] = y + u2Row[x]; - } + for (int x = 0; x < u1.cols; ++x) + { + map1Row[x] = x + u1Row[x]; + map2Row[x] = y + u2Row[x]; } } +} - void buildFlowMap(const Mat_& u1, const Mat_& u2, Mat_& map1, Mat_& map2) - { - CV_DbgAssert( u2.size() == u1.size() ); - CV_DbgAssert( map1.size() == u1.size() ); - CV_DbgAssert( map2.size() == u1.size() ); +void buildFlowMap(const Mat_& u1, const Mat_& u2, Mat_& map1, Mat_& map2) +{ + CV_DbgAssert( u2.size() == u1.size() ); + CV_DbgAssert( map1.size() == u1.size() ); + CV_DbgAssert( map2.size() == u1.size() ); - BuildFlowMapBody body; + BuildFlowMapBody body; - body.u1 = u1; - body.u2 = u2; - body.map1 = map1; - body.map2 = map2; + body.u1 = u1; + body.u2 = u2; + body.map1 = map1; + body.map2 = map2; - parallel_for_(Range(0, u1.rows), body); - } + parallel_for_(Range(0, u1.rows), body); +} - //////////////////////////////////////////////////////////// - // centeredGradient +//////////////////////////////////////////////////////////// +// centeredGradient - struct CenteredGradientBody : ParallelLoopBody - { - void operator() (const Range& range) const; +struct CenteredGradientBody : ParallelLoopBody +{ + void operator() (const Range& range) const; - Mat_ src; - mutable Mat_ dx; - mutable Mat_ dy; - }; + Mat_ src; + mutable Mat_ dx; + mutable Mat_ dy; +}; - void CenteredGradientBody::operator() (const Range& range) const - { - const int last_col = src.cols - 1; +void CenteredGradientBody::operator() (const Range& range) const +{ + const int last_col = src.cols - 1; - for (int y = range.start; y < range.end; ++y) - { - const float* srcPrevRow = src[y - 1]; - const float* srcCurRow = src[y]; - const float* srcNextRow = src[y + 1]; + for (int y = range.start; y < range.end; ++y) + { + const float* srcPrevRow = src[y - 1]; + const float* srcCurRow = src[y]; + const float* srcNextRow = src[y + 1]; - float* dxRow = dx[y]; - float* dyRow = dy[y]; + float* dxRow = dx[y]; + float* dyRow = dy[y]; - for (int x = 1; x < last_col; ++x) - { - dxRow[x] = 0.5f * (srcCurRow[x + 1] - srcCurRow[x - 1]); - dyRow[x] = 0.5f * (srcNextRow[x] - srcPrevRow[x]); - } + for (int x = 1; x < last_col; ++x) + { + dxRow[x] = 0.5f * (srcCurRow[x + 1] - srcCurRow[x - 1]); + dyRow[x] = 0.5f * (srcNextRow[x] - srcPrevRow[x]); } } +} - void centeredGradient(const Mat_& src, Mat_& dx, Mat_& dy) - { - CV_DbgAssert( src.rows > 2 && src.cols > 2 ); - CV_DbgAssert( dx.size() == src.size() ); - CV_DbgAssert( dy.size() == src.size() ); +void centeredGradient(const Mat_& src, Mat_& dx, Mat_& dy) +{ + CV_DbgAssert( src.rows > 2 && src.cols > 2 ); + CV_DbgAssert( dx.size() == src.size() ); + CV_DbgAssert( dy.size() == src.size() ); - const int last_row = src.rows - 1; - const int last_col = src.cols - 1; + const int last_row = src.rows - 1; + const int last_col = src.cols - 1; - // compute the gradient on the center body of the image - { - CenteredGradientBody body; + // compute the gradient on the center body of the image + { + CenteredGradientBody body; - body.src = src; - body.dx = dx; - body.dy = dy; + body.src = src; + body.dx = dx; + body.dy = dy; - parallel_for_(Range(1, last_row), body); - } + parallel_for_(Range(1, last_row), body); + } - // compute the gradient on the first and last rows - for (int x = 1; x < last_col; ++x) - { - dx(0, x) = 0.5f * (src(0, x + 1) - src(0, x - 1)); - dy(0, x) = 0.5f * (src(1, x) - src(0, x)); + // compute the gradient on the first and last rows + for (int x = 1; x < last_col; ++x) + { + dx(0, x) = 0.5f * (src(0, x + 1) - src(0, x - 1)); + dy(0, x) = 0.5f * (src(1, x) - src(0, x)); - dx(last_row, x) = 0.5f * (src(last_row, x + 1) - src(last_row, x - 1)); - dy(last_row, x) = 0.5f * (src(last_row, x) - src(last_row - 1, x)); - } + dx(last_row, x) = 0.5f * (src(last_row, x + 1) - src(last_row, x - 1)); + dy(last_row, x) = 0.5f * (src(last_row, x) - src(last_row - 1, x)); + } - // compute the gradient on the first and last columns - for (int y = 1; y < last_row; ++y) - { - dx(y, 0) = 0.5f * (src(y, 1) - src(y, 0)); - dy(y, 0) = 0.5f * (src(y + 1, 0) - src(y - 1, 0)); + // compute the gradient on the first and last columns + for (int y = 1; y < last_row; ++y) + { + dx(y, 0) = 0.5f * (src(y, 1) - src(y, 0)); + dy(y, 0) = 0.5f * (src(y + 1, 0) - src(y - 1, 0)); - dx(y, last_col) = 0.5f * (src(y, last_col) - src(y, last_col - 1)); - dy(y, last_col) = 0.5f * (src(y + 1, last_col) - src(y - 1, last_col)); - } + dx(y, last_col) = 0.5f * (src(y, last_col) - src(y, last_col - 1)); + dy(y, last_col) = 0.5f * (src(y + 1, last_col) - src(y - 1, last_col)); + } - // compute the gradient at the four corners - dx(0, 0) = 0.5f * (src(0, 1) - src(0, 0)); - dy(0, 0) = 0.5f * (src(1, 0) - src(0, 0)); + // compute the gradient at the four corners + dx(0, 0) = 0.5f * (src(0, 1) - src(0, 0)); + dy(0, 0) = 0.5f * (src(1, 0) - src(0, 0)); - dx(0, last_col) = 0.5f * (src(0, last_col) - src(0, last_col - 1)); - dy(0, last_col) = 0.5f * (src(1, last_col) - src(0, last_col)); + dx(0, last_col) = 0.5f * (src(0, last_col) - src(0, last_col - 1)); + dy(0, last_col) = 0.5f * (src(1, last_col) - src(0, last_col)); - dx(last_row, 0) = 0.5f * (src(last_row, 1) - src(last_row, 0)); - dy(last_row, 0) = 0.5f * (src(last_row, 0) - src(last_row - 1, 0)); + dx(last_row, 0) = 0.5f * (src(last_row, 1) - src(last_row, 0)); + dy(last_row, 0) = 0.5f * (src(last_row, 0) - src(last_row - 1, 0)); - dx(last_row, last_col) = 0.5f * (src(last_row, last_col) - src(last_row, last_col - 1)); - dy(last_row, last_col) = 0.5f * (src(last_row, last_col) - src(last_row - 1, last_col)); - } + dx(last_row, last_col) = 0.5f * (src(last_row, last_col) - src(last_row, last_col - 1)); + dy(last_row, last_col) = 0.5f * (src(last_row, last_col) - src(last_row - 1, last_col)); +} - //////////////////////////////////////////////////////////// - // forwardGradient +//////////////////////////////////////////////////////////// +// forwardGradient - struct ForwardGradientBody : ParallelLoopBody - { - void operator() (const Range& range) const; +struct ForwardGradientBody : ParallelLoopBody +{ + void operator() (const Range& range) const; - Mat_ src; - mutable Mat_ dx; - mutable Mat_ dy; - }; + Mat_ src; + mutable Mat_ dx; + mutable Mat_ dy; +}; - void ForwardGradientBody::operator() (const Range& range) const - { - const int last_col = src.cols - 1; +void ForwardGradientBody::operator() (const Range& range) const +{ + const int last_col = src.cols - 1; - for (int y = range.start; y < range.end; ++y) - { - const float* srcCurRow = src[y]; - const float* srcNextRow = src[y + 1]; + for (int y = range.start; y < range.end; ++y) + { + const float* srcCurRow = src[y]; + const float* srcNextRow = src[y + 1]; - float* dxRow = dx[y]; - float* dyRow = dy[y]; + float* dxRow = dx[y]; + float* dyRow = dy[y]; - for (int x = 0; x < last_col; ++x) - { - dxRow[x] = srcCurRow[x + 1] - srcCurRow[x]; - dyRow[x] = srcNextRow[x] - srcCurRow[x]; - } + for (int x = 0; x < last_col; ++x) + { + dxRow[x] = srcCurRow[x + 1] - srcCurRow[x]; + dyRow[x] = srcNextRow[x] - srcCurRow[x]; } } +} - void forwardGradient(const Mat_& src, Mat_& dx, Mat_& dy) - { - CV_DbgAssert( src.rows > 2 && src.cols > 2 ); - CV_DbgAssert( dx.size() == src.size() ); - CV_DbgAssert( dy.size() == src.size() ); - - const int last_row = src.rows - 1; - const int last_col = src.cols - 1; +void forwardGradient(const Mat_& src, Mat_& dx, Mat_& dy) +{ + CV_DbgAssert( src.rows > 2 && src.cols > 2 ); + CV_DbgAssert( dx.size() == src.size() ); + CV_DbgAssert( dy.size() == src.size() ); - // compute the gradient on the central body of the image - { - ForwardGradientBody body; + const int last_row = src.rows - 1; + const int last_col = src.cols - 1; - body.src = src; - body.dx = dx; - body.dy = dy; + // compute the gradient on the central body of the image + { + ForwardGradientBody body; - parallel_for_(Range(0, last_row), body); - } + body.src = src; + body.dx = dx; + body.dy = dy; - // compute the gradient on the last row - for (int x = 0; x < last_col; ++x) - { - dx(last_row, x) = src(last_row, x + 1) - src(last_row, x); - dy(last_row, x) = 0.0f; - } + parallel_for_(Range(0, last_row), body); + } - // compute the gradient on the last column - for (int y = 0; y < last_row; ++y) - { - dx(y, last_col) = 0.0f; - dy(y, last_col) = src(y + 1, last_col) - src(y, last_col); - } + // compute the gradient on the last row + for (int x = 0; x < last_col; ++x) + { + dx(last_row, x) = src(last_row, x + 1) - src(last_row, x); + dy(last_row, x) = 0.0f; + } - dx(last_row, last_col) = 0.0f; - dy(last_row, last_col) = 0.0f; + // compute the gradient on the last column + for (int y = 0; y < last_row; ++y) + { + dx(y, last_col) = 0.0f; + dy(y, last_col) = src(y + 1, last_col) - src(y, last_col); } - //////////////////////////////////////////////////////////// - // divergence + dx(last_row, last_col) = 0.0f; + dy(last_row, last_col) = 0.0f; +} + +//////////////////////////////////////////////////////////// +// divergence - struct DivergenceBody : ParallelLoopBody - { - void operator() (const Range& range) const; +struct DivergenceBody : ParallelLoopBody +{ + void operator() (const Range& range) const; - Mat_ v1; - Mat_ v2; - mutable Mat_ div; - }; + Mat_ v1; + Mat_ v2; + mutable Mat_ div; +}; - void DivergenceBody::operator() (const Range& range) const +void DivergenceBody::operator() (const Range& range) const +{ + for (int y = range.start; y < range.end; ++y) { - for (int y = range.start; y < range.end; ++y) - { - const float* v1Row = v1[y]; - const float* v2PrevRow = v2[y - 1]; - const float* v2CurRow = v2[y]; + const float* v1Row = v1[y]; + const float* v2PrevRow = v2[y - 1]; + const float* v2CurRow = v2[y]; - float* divRow = div[y]; + float* divRow = div[y]; - for(int x = 1; x < v1.cols; ++x) - { - const float v1x = v1Row[x] - v1Row[x - 1]; - const float v2y = v2CurRow[x] - v2PrevRow[x]; + for(int x = 1; x < v1.cols; ++x) + { + const float v1x = v1Row[x] - v1Row[x - 1]; + const float v2y = v2CurRow[x] - v2PrevRow[x]; - divRow[x] = v1x + v2y; - } + divRow[x] = v1x + v2y; } } +} + +void divergence(const Mat_& v1, const Mat_& v2, Mat_& div) +{ + CV_DbgAssert( v1.rows > 2 && v1.cols > 2 ); + CV_DbgAssert( v2.size() == v1.size() ); + CV_DbgAssert( div.size() == v1.size() ); - void divergence(const Mat_& v1, const Mat_& v2, Mat_& div) { - CV_DbgAssert( v1.rows > 2 && v1.cols > 2 ); - CV_DbgAssert( v2.size() == v1.size() ); - CV_DbgAssert( div.size() == v1.size() ); + DivergenceBody body; - { - DivergenceBody body; + body.v1 = v1; + body.v2 = v2; + body.div = div; - body.v1 = v1; - body.v2 = v2; - body.div = div; + parallel_for_(Range(1, v1.rows), body); + } - parallel_for_(Range(1, v1.rows), body); - } + // compute the divergence on the first row + for(int x = 1; x < v1.cols; ++x) + div(0, x) = v1(0, x) - v1(0, x - 1) + v2(0, x); - // compute the divergence on the first row - for(int x = 1; x < v1.cols; ++x) - div(0, x) = v1(0, x) - v1(0, x - 1) + v2(0, x); + // compute the divergence on the first column + for (int y = 1; y < v1.rows; ++y) + div(y, 0) = v1(y, 0) + v2(y, 0) - v2(y - 1, 0); - // compute the divergence on the first column - for (int y = 1; y < v1.rows; ++y) - div(y, 0) = v1(y, 0) + v2(y, 0) - v2(y - 1, 0); - - div(0, 0) = v1(0, 0) + v2(0, 0); - } + div(0, 0) = v1(0, 0) + v2(0, 0); +} - //////////////////////////////////////////////////////////// - // calcGradRho +//////////////////////////////////////////////////////////// +// calcGradRho - struct CalcGradRhoBody : ParallelLoopBody - { - void operator() (const Range& range) const; - - Mat_ I0; - Mat_ I1w; - Mat_ I1wx; - Mat_ I1wy; - Mat_ u1; - Mat_ u2; - mutable Mat_ grad; - mutable Mat_ rho_c; - }; - - void CalcGradRhoBody::operator() (const Range& range) const +struct CalcGradRhoBody : ParallelLoopBody +{ + void operator() (const Range& range) const; + + Mat_ I0; + Mat_ I1w; + Mat_ I1wx; + Mat_ I1wy; + Mat_ u1; + Mat_ u2; + mutable Mat_ grad; + mutable Mat_ rho_c; +}; + +void CalcGradRhoBody::operator() (const Range& range) const +{ + for (int y = range.start; y < range.end; ++y) { - for (int y = range.start; y < range.end; ++y) - { - const float* I0Row = I0[y]; - const float* I1wRow = I1w[y]; - const float* I1wxRow = I1wx[y]; - const float* I1wyRow = I1wy[y]; - const float* u1Row = u1[y]; - const float* u2Row = u2[y]; + const float* I0Row = I0[y]; + const float* I1wRow = I1w[y]; + const float* I1wxRow = I1wx[y]; + const float* I1wyRow = I1wy[y]; + const float* u1Row = u1[y]; + const float* u2Row = u2[y]; - float* gradRow = grad[y]; - float* rhoRow = rho_c[y]; + float* gradRow = grad[y]; + float* rhoRow = rho_c[y]; - for (int x = 0; x < I0.cols; ++x) - { - const float Ix2 = I1wxRow[x] * I1wxRow[x]; - const float Iy2 = I1wyRow[x] * I1wyRow[x]; + for (int x = 0; x < I0.cols; ++x) + { + const float Ix2 = I1wxRow[x] * I1wxRow[x]; + const float Iy2 = I1wyRow[x] * I1wyRow[x]; - // store the |Grad(I1)|^2 - gradRow[x] = Ix2 + Iy2; + // store the |Grad(I1)|^2 + gradRow[x] = Ix2 + Iy2; - // compute the constant part of the rho function - rhoRow[x] = (I1wRow[x] - I1wxRow[x] * u1Row[x] - I1wyRow[x] * u2Row[x] - I0Row[x]); - } + // compute the constant part of the rho function + rhoRow[x] = (I1wRow[x] - I1wxRow[x] * u1Row[x] - I1wyRow[x] * u2Row[x] - I0Row[x]); } } +} - void calcGradRho(const Mat_& I0, const Mat_& I1w, const Mat_& I1wx, const Mat_& I1wy, const Mat_& u1, const Mat_& u2, - Mat_& grad, Mat_& rho_c) - { - CV_DbgAssert( I1w.size() == I0.size() ); - CV_DbgAssert( I1wx.size() == I0.size() ); - CV_DbgAssert( I1wy.size() == I0.size() ); - CV_DbgAssert( u1.size() == I0.size() ); - CV_DbgAssert( u2.size() == I0.size() ); - CV_DbgAssert( grad.size() == I0.size() ); - CV_DbgAssert( rho_c.size() == I0.size() ); - - CalcGradRhoBody body; - - body.I0 = I0; - body.I1w = I1w; - body.I1wx = I1wx; - body.I1wy = I1wy; - body.u1 = u1; - body.u2 = u2; - body.grad = grad; - body.rho_c = rho_c; - - parallel_for_(Range(0, I0.rows), body); - } +void calcGradRho(const Mat_& I0, const Mat_& I1w, const Mat_& I1wx, const Mat_& I1wy, const Mat_& u1, const Mat_& u2, + Mat_& grad, Mat_& rho_c) +{ + CV_DbgAssert( I1w.size() == I0.size() ); + CV_DbgAssert( I1wx.size() == I0.size() ); + CV_DbgAssert( I1wy.size() == I0.size() ); + CV_DbgAssert( u1.size() == I0.size() ); + CV_DbgAssert( u2.size() == I0.size() ); + CV_DbgAssert( grad.size() == I0.size() ); + CV_DbgAssert( rho_c.size() == I0.size() ); + + CalcGradRhoBody body; + + body.I0 = I0; + body.I1w = I1w; + body.I1wx = I1wx; + body.I1wy = I1wy; + body.u1 = u1; + body.u2 = u2; + body.grad = grad; + body.rho_c = rho_c; + + parallel_for_(Range(0, I0.rows), body); +} - //////////////////////////////////////////////////////////// - // estimateV +//////////////////////////////////////////////////////////// +// estimateV - struct EstimateVBody : ParallelLoopBody - { - void operator() (const Range& range) const; - - Mat_ I1wx; - Mat_ I1wy; - Mat_ u1; - Mat_ u2; - Mat_ grad; - Mat_ rho_c; - mutable Mat_ v1; - mutable Mat_ v2; - float l_t; - }; - - void EstimateVBody::operator() (const Range& range) const +struct EstimateVBody : ParallelLoopBody +{ + void operator() (const Range& range) const; + + Mat_ I1wx; + Mat_ I1wy; + Mat_ u1; + Mat_ u2; + Mat_ grad; + Mat_ rho_c; + mutable Mat_ v1; + mutable Mat_ v2; + float l_t; +}; + +void EstimateVBody::operator() (const Range& range) const +{ + for (int y = range.start; y < range.end; ++y) { - 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* gradRow = grad[y]; + const float* rhoRow = rho_c[y]; + + float* v1Row = v1[y]; + float* v2Row = v2[y]; + + for (int x = 0; x < I1wx.cols; ++x) { - const float* I1wxRow = I1wx[y]; - const float* I1wyRow = I1wy[y]; - const float* u1Row = u1[y]; - const float* u2Row = u2[y]; - const float* gradRow = grad[y]; - const float* rhoRow = rho_c[y]; + const float rho = rhoRow[x] + (I1wxRow[x] * u1Row[x] + I1wyRow[x] * u2Row[x]); - float* v1Row = v1[y]; - float* v2Row = v2[y]; + float d1 = 0.0f; + float d2 = 0.0f; - for (int x = 0; x < I1wx.cols; ++x) + if (rho < -l_t * gradRow[x]) { - const float rho = rhoRow[x] + (I1wxRow[x] * u1Row[x] + I1wyRow[x] * u2Row[x]); - - float d1 = 0.0f; - float d2 = 0.0f; - - if (rho < -l_t * gradRow[x]) - { - d1 = l_t * I1wxRow[x]; - d2 = l_t * I1wyRow[x]; - } - else if (rho > l_t * gradRow[x]) - { - d1 = -l_t * I1wxRow[x]; - d2 = -l_t * I1wyRow[x]; - } - else if (gradRow[x] > numeric_limits::epsilon()) - { - float fi = -rho / gradRow[x]; - d1 = fi * I1wxRow[x]; - d2 = fi * I1wyRow[x]; - } - - v1Row[x] = u1Row[x] + d1; - v2Row[x] = u2Row[x] + d2; + d1 = l_t * I1wxRow[x]; + d2 = l_t * I1wyRow[x]; } + else if (rho > l_t * gradRow[x]) + { + d1 = -l_t * I1wxRow[x]; + d2 = -l_t * I1wyRow[x]; + } + else if (gradRow[x] > numeric_limits::epsilon()) + { + float fi = -rho / gradRow[x]; + d1 = fi * I1wxRow[x]; + d2 = fi * I1wyRow[x]; + } + + v1Row[x] = u1Row[x] + d1; + v2Row[x] = u2Row[x] + d2; } } +} - void estimateV(const Mat_& I1wx, const Mat_& I1wy, const Mat_& u1, const Mat_& u2, const Mat_& grad, const Mat_& rho_c, - Mat_& v1, Mat_& v2, float l_t) - { - CV_DbgAssert( I1wy.size() == I1wx.size() ); - CV_DbgAssert( u1.size() == I1wx.size() ); - CV_DbgAssert( u2.size() == I1wx.size() ); - CV_DbgAssert( grad.size() == I1wx.size() ); - CV_DbgAssert( rho_c.size() == I1wx.size() ); - CV_DbgAssert( v1.size() == I1wx.size() ); - CV_DbgAssert( v2.size() == I1wx.size() ); - - EstimateVBody body; - - body.I1wx = I1wx; - body.I1wy = I1wy; - body.u1 = u1; - body.u2 = u2; - body.grad = grad; - body.rho_c = rho_c; - body.v1 = v1; - body.v2 = v2; - body.l_t = l_t; - - parallel_for_(Range(0, I1wx.rows), body); - } +void estimateV(const Mat_& I1wx, const Mat_& I1wy, const Mat_& u1, const Mat_& u2, const Mat_& grad, const Mat_& rho_c, + Mat_& v1, Mat_& v2, float l_t) +{ + CV_DbgAssert( I1wy.size() == I1wx.size() ); + CV_DbgAssert( u1.size() == I1wx.size() ); + CV_DbgAssert( u2.size() == I1wx.size() ); + CV_DbgAssert( grad.size() == I1wx.size() ); + CV_DbgAssert( rho_c.size() == I1wx.size() ); + CV_DbgAssert( v1.size() == I1wx.size() ); + CV_DbgAssert( v2.size() == I1wx.size() ); + + EstimateVBody body; + + body.I1wx = I1wx; + body.I1wy = I1wy; + body.u1 = u1; + body.u2 = u2; + body.grad = grad; + body.rho_c = rho_c; + body.v1 = v1; + body.v2 = v2; + body.l_t = l_t; + + parallel_for_(Range(0, I1wx.rows), body); +} - //////////////////////////////////////////////////////////// - // estimateU +//////////////////////////////////////////////////////////// +// estimateU - float estimateU(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_& div_p1, const Mat_& div_p2, Mat_& u1, Mat_& u2, float theta) +{ + CV_DbgAssert( v2.size() == v1.size() ); + CV_DbgAssert( div_p1.size() == v1.size() ); + CV_DbgAssert( div_p2.size() == v1.size() ); + CV_DbgAssert( u1.size() == v1.size() ); + CV_DbgAssert( u2.size() == v1.size() ); + + float error = 0.0f; + for (int y = 0; y < v1.rows; ++y) { - CV_DbgAssert( v2.size() == v1.size() ); - CV_DbgAssert( div_p1.size() == v1.size() ); - CV_DbgAssert( div_p2.size() == v1.size() ); - CV_DbgAssert( u1.size() == v1.size() ); - CV_DbgAssert( u2.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* divP1Row = div_p1[y]; - const float* divP2Row = div_p2[y]; + const float* v1Row = v1[y]; + const float* v2Row = v2[y]; + const float* divP1Row = div_p1[y]; + const float* divP2Row = div_p2[y]; - float* u1Row = u1[y]; - float* u2Row = u2[y]; + float* u1Row = u1[y]; + float* u2Row = u2[y]; - for (int x = 0; x < v1.cols; ++x) - { - const float u1k = u1Row[x]; - const float u2k = u2Row[x]; + for (int x = 0; x < v1.cols; ++x) + { + const float u1k = u1Row[x]; + const float u2k = u2Row[x]; - u1Row[x] = v1Row[x] + theta * divP1Row[x]; - u2Row[x] = v2Row[x] + theta * divP2Row[x]; + 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); - } + error += (u1Row[x] - u1k) * (u1Row[x] - u1k) + (u2Row[x] - u2k) * (u2Row[x] - u2k); } - - return error; } - //////////////////////////////////////////////////////////// - // estimateDualVariables + return error; +} + +//////////////////////////////////////////////////////////// +// estimateDualVariables - struct EstimateDualVariablesBody : ParallelLoopBody - { - void operator() (const Range& range) const; - - Mat_ u1x; - Mat_ u1y; - Mat_ u2x; - Mat_ u2y; - mutable Mat_ p11; - mutable Mat_ p12; - mutable Mat_ p21; - mutable Mat_ p22; - float taut; - }; - - void EstimateDualVariablesBody::operator() (const Range& range) const +struct EstimateDualVariablesBody : ParallelLoopBody +{ + void operator() (const Range& range) const; + + Mat_ u1x; + Mat_ u1y; + Mat_ u2x; + Mat_ u2y; + mutable Mat_ p11; + mutable Mat_ p12; + mutable Mat_ p21; + mutable Mat_ p22; + float taut; +}; + +void EstimateDualVariablesBody::operator() (const Range& range) const +{ + for (int y = range.start; y < range.end; ++y) { - for (int y = range.start; y < range.end; ++y) - { - const float* u1xRow = u1x[y]; - const float* u1yRow = u1y[y]; - const float* u2xRow = u2x[y]; - const float* u2yRow = u2y[y]; + const float* u1xRow = u1x[y]; + const float* u1yRow = u1y[y]; + const float* u2xRow = u2x[y]; + const float* u2yRow = u2y[y]; - float* p11Row = p11[y]; - float* p12Row = p12[y]; - float* p21Row = p21[y]; - float* p22Row = p22[y]; + float* p11Row = p11[y]; + float* p12Row = p12[y]; + float* p21Row = p21[y]; + float* p22Row = p22[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])); + 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 ng1 = 1.0f + taut * g1; - const float ng2 = 1.0f + taut * g2; + const float ng1 = 1.0f + taut * g1; + const float ng2 = 1.0f + taut * g2; - 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; - } + 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; } } +} - void estimateDualVariables(const Mat_& u1x, const Mat_& u1y, const Mat_& u2x, const Mat_& u2y, - Mat_& p11, Mat_& p12, Mat_& p21, Mat_& p22, float taut) - { - CV_DbgAssert( u1y.size() == u1x.size() ); - CV_DbgAssert( u2x.size() == u1x.size() ); - CV_DbgAssert( u2y.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() ); - - EstimateDualVariablesBody body; - - body.u1x = u1x; - body.u1y = u1y; - body.u2x = u2x; - body.u2y = u2y; - body.p11 = p11; - body.p12 = p12; - body.p21 = p21; - body.p22 = p22; - body.taut = taut; - - parallel_for_(Range(0, u1x.rows), body); - } +void estimateDualVariables(const Mat_& u1x, const Mat_& u1y, const Mat_& u2x, const Mat_& u2y, + Mat_& p11, Mat_& p12, Mat_& p21, Mat_& p22, float taut) +{ + CV_DbgAssert( u1y.size() == u1x.size() ); + CV_DbgAssert( u2x.size() == u1x.size() ); + CV_DbgAssert( u2y.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() ); + + EstimateDualVariablesBody body; + + body.u1x = u1x; + body.u1y = u1y; + body.u2x = u2x; + body.u2y = u2y; + body.p11 = p11; + body.p12 = p12; + body.p21 = p21; + body.p22 = p22; + body.taut = taut; + + parallel_for_(Range(0, u1x.rows), body); } -void cv::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) { const float scaledEpsilon = static_cast(epsilon * epsilon * I0.size().area()); @@ -818,21 +875,12 @@ void cv::OpticalFlowDual_TVL1::procOneScale(const Mat_& I0, const Mat_ void releaseVector(vector& v) - { - vector empty; - empty.swap(v); - } -} - -void cv::OpticalFlowDual_TVL1::collectGarbage() +void OpticalFlowDual_TVL1::collectGarbage() { - releaseVector(I0s); - releaseVector(I1s); - releaseVector(u1s); - releaseVector(u2s); + I0s.clear(); + I1s.clear(); + u1s.clear(); + u2s.clear(); I1x_buf.release(); I1y_buf.release(); @@ -863,3 +911,27 @@ void cv::OpticalFlowDual_TVL1::collectGarbage() u2x_buf.release(); u2y_buf.release(); } + +CV_INIT_ALGORITHM(OpticalFlowDual_TVL1, "DenseOpticalFlow.DualTVL1", + obj.info()->addParam(obj, "tau", obj.tau, false, 0, 0, + "Time step of the numerical scheme"); + obj.info()->addParam(obj, "lambda", obj.lambda, false, 0, 0, + "Weight parameter for the data term, attachment parameter"); + obj.info()->addParam(obj, "theta", obj.theta, false, 0, 0, + "Weight parameter for (u - v)^2, tightness parameter"); + obj.info()->addParam(obj, "nscales", obj.nscales, false, 0, 0, + "Number of scales used to create the pyramid of images"); + obj.info()->addParam(obj, "warps", obj.warps, false, 0, 0, + "Number of warpings per scale"); + obj.info()->addParam(obj, "epsilon", obj.epsilon, false, 0, 0, + "Stopping criterion threshold used in the numerical scheme, which is a trade-off between precision and running time"); + obj.info()->addParam(obj, "iterations", obj.iterations, false, 0, 0, + "Stopping criterion iterations number used in the numerical scheme"); + obj.info()->addParam(obj, "useInitialFlow", obj.useInitialFlow)); + +} // namespace + +Ptr cv::createOptFlow_DualTVL1() +{ + return new OpticalFlowDual_TVL1; +} diff --git a/modules/video/test/test_tvl1optflow.cpp b/modules/video/test/test_tvl1optflow.cpp index b5688d3..804eae8 100644 --- a/modules/video/test/test_tvl1optflow.cpp +++ b/modules/video/test/test_tvl1optflow.cpp @@ -152,9 +152,9 @@ TEST(Video_calcOpticalFlowDual_TVL1, Regression) ASSERT_FALSE(frame2.empty()); Mat_ flow; - OpticalFlowDual_TVL1 tvl1; + Ptr tvl1 = createOptFlow_DualTVL1(); - tvl1(frame1, frame2, flow); + tvl1->calc(frame1, frame2, flow); #ifdef DUMP writeOpticalFlowToFile(flow, gold_flow_path); diff --git a/samples/cpp/tvl1_optical_flow.cpp b/samples/cpp/tvl1_optical_flow.cpp index d9a57a2..4f9a02b 100644 --- a/samples/cpp/tvl1_optical_flow.cpp +++ b/samples/cpp/tvl1_optical_flow.cpp @@ -173,10 +173,10 @@ int main(int argc, const char* argv[]) } Mat_ flow; - OpticalFlowDual_TVL1 tvl1; + Ptr tvl1 = createOptFlow_DualTVL1(); const double start = (double)getTickCount(); - tvl1(frame0, frame1, flow); + tvl1->calc(frame0, frame1, flow); const double timeSec = (getTickCount() - start) / getTickFrequency(); cout << "calcOpticalFlowDual_TVL1 : " << timeSec << " sec" << endl; -- 2.7.4