Commit a3a09cf4 authored by Vladislav Vinogradov's avatar Vladislav Vinogradov

refactored OpticalFlowDual_TVL1:

* added DenseOpticalFlow interface
* moved OpticalFlowDual_TVL1 to src folder
parent 2181a41a
...@@ -431,13 +431,13 @@ PERF_TEST_P(ImagePair, Video_OpticalFlowDual_TVL1, ...@@ -431,13 +431,13 @@ PERF_TEST_P(ImagePair, Video_OpticalFlowDual_TVL1,
{ {
cv::Mat flow; cv::Mat flow;
cv::OpticalFlowDual_TVL1 alg; cv::Ptr<cv::DenseOpticalFlow> alg = cv::createOptFlow_DualTVL1();
alg(frame0, frame1, flow); alg->calc(frame0, frame1, flow);
TEST_CYCLE() TEST_CYCLE()
{ {
alg(frame0, frame1, flow); alg->calc(frame0, frame1, flow);
} }
CPU_SANITY_CHECK(flow); CPU_SANITY_CHECK(flow);
......
...@@ -431,9 +431,9 @@ GPU_TEST_P(OpticalFlowDual_TVL1, Accuracy) ...@@ -431,9 +431,9 @@ GPU_TEST_P(OpticalFlowDual_TVL1, Accuracy)
cv::gpu::GpuMat d_flowy = createMat(frame0.size(), CV_32FC1, useRoi); cv::gpu::GpuMat d_flowy = createMat(frame0.size(), CV_32FC1, useRoi);
d_alg(loadMat(frame0, useRoi), loadMat(frame1, useRoi), d_flowx, d_flowy); d_alg(loadMat(frame0, useRoi), loadMat(frame1, useRoi), d_flowx, d_flowy);
cv::OpticalFlowDual_TVL1 alg; cv::Ptr<cv::DenseOpticalFlow> alg = cv::createOptFlow_DualTVL1();
cv::Mat flow; cv::Mat flow;
alg(frame0, frame1, flow); alg->calc(frame0, frame1, flow);
cv::Mat gold[2]; cv::Mat gold[2];
cv::split(flow, gold); cv::split(flow, gold);
......
...@@ -643,11 +643,11 @@ See [Tao2012]_. And site of project - http://graphics.berkeley.edu/papers/Tao-SA ...@@ -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. "Dual TV L1" Optical Flow Algorithm.
.. ocv:class:: OpticalFlowDual_TVL12 .. ocv:function:: Ptr<DenseOpticalFlow> createOptFlow_DualTVL1()
The class implements the "Dual TV L1" optical flow algorithm described in [Zach2007]_ and [Javier2012]_ . 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 ...@@ -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. 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. :param prev: first 8-bit single-channel input image.
...@@ -699,11 +699,11 @@ Calculates an optical flow. ...@@ -699,11 +699,11 @@ Calculates an optical flow.
OpticalFlowDual_TVL1::collectGarbage DenseOpticalFlow::collectGarbage
------------------------------------ --------------------------------
Releases all inner buffers. Releases all inner buffers.
.. ocv:function:: void OpticalFlowDual_TVL1::collectGarbage() .. ocv:function:: void DenseOpticalFlow::collectGarbage()
......
...@@ -352,104 +352,19 @@ CV_EXPORTS_W void calcOpticalFlowSF(Mat& from, ...@@ -352,104 +352,19 @@ CV_EXPORTS_W void calcOpticalFlowSF(Mat& from,
double upscale_sigma_color, double upscale_sigma_color,
double speed_up_thr); 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 // Implementation of the Zach, Pock and Bischof Dual TV-L1 Optical Flow method
// //
// see reference: // see reference:
// [1] C. Zach, T. Pock and H. Bischof, "A Duality Based Approach for Realtime TV-L1 Optical Flow". // [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". // [2] Javier Sanchez, Enric Meinhardt-Llopis and Gabriele Facciolo. "TV-L1 Optical Flow Estimation".
class CV_EXPORTS OpticalFlowDual_TVL1 CV_EXPORTS Ptr<DenseOpticalFlow> createOptFlow_DualTVL1();
{
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_<float>& I0, const Mat_<float>& I1, Mat_<float>& u1, Mat_<float>& u2);
std::vector<Mat_<float> > I0s;
std::vector<Mat_<float> > I1s;
std::vector<Mat_<float> > u1s;
std::vector<Mat_<float> > u2s;
Mat_<float> I1x_buf;
Mat_<float> I1y_buf;
Mat_<float> flowMap1_buf;
Mat_<float> flowMap2_buf;
Mat_<float> I1w_buf;
Mat_<float> I1wx_buf;
Mat_<float> I1wy_buf;
Mat_<float> grad_buf;
Mat_<float> rho_c_buf;
Mat_<float> v1_buf;
Mat_<float> v2_buf;
Mat_<float> p11_buf;
Mat_<float> p12_buf;
Mat_<float> p21_buf;
Mat_<float> p22_buf;
Mat_<float> div_p1_buf;
Mat_<float> div_p2_buf;
Mat_<float> u1x_buf;
Mat_<float> u1y_buf;
Mat_<float> u2x_buf;
Mat_<float> u2y_buf;
};
} }
......
...@@ -22,12 +22,9 @@ PERF_TEST_P(ImagePair, OpticalFlowDual_TVL1, testing::Values(impair("cv/optflow/ ...@@ -22,12 +22,9 @@ PERF_TEST_P(ImagePair, OpticalFlowDual_TVL1, testing::Values(impair("cv/optflow/
Mat flow; Mat flow;
OpticalFlowDual_TVL1 tvl1; Ptr<DenseOpticalFlow> tvl1 = createOptFlow_DualTVL1();
TEST_CYCLE() TEST_CYCLE_N(10) tvl1->calc(frame1, frame2, flow);
{
tvl1(frame1, frame2, flow);
}
SANITY_CHECK(flow, 0.5); SANITY_CHECK(flow, 0.5);
} }
...@@ -77,7 +77,67 @@ ...@@ -77,7 +77,67 @@
using namespace std; using namespace std;
using namespace cv; 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_<float>& I0, const Mat_<float>& I1, Mat_<float>& u1, Mat_<float>& u2);
std::vector<Mat_<float> > I0s;
std::vector<Mat_<float> > I1s;
std::vector<Mat_<float> > u1s;
std::vector<Mat_<float> > u2s;
Mat_<float> I1x_buf;
Mat_<float> I1y_buf;
Mat_<float> flowMap1_buf;
Mat_<float> flowMap2_buf;
Mat_<float> I1w_buf;
Mat_<float> I1wx_buf;
Mat_<float> I1wy_buf;
Mat_<float> grad_buf;
Mat_<float> rho_c_buf;
Mat_<float> v1_buf;
Mat_<float> v2_buf;
Mat_<float> p11_buf;
Mat_<float> p12_buf;
Mat_<float> p21_buf;
Mat_<float> p22_buf;
Mat_<float> div_p1_buf;
Mat_<float> div_p2_buf;
Mat_<float> u1x_buf;
Mat_<float> u1y_buf;
Mat_<float> u2x_buf;
Mat_<float> u2y_buf;
};
OpticalFlowDual_TVL1::OpticalFlowDual_TVL1()
{ {
tau = 0.25; tau = 0.25;
lambda = 0.15; lambda = 0.15;
...@@ -89,7 +149,7 @@ cv::OpticalFlowDual_TVL1::OpticalFlowDual_TVL1() ...@@ -89,7 +149,7 @@ cv::OpticalFlowDual_TVL1::OpticalFlowDual_TVL1()
useInitialFlow = false; 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 I0 = _I0.getMat();
Mat I1 = _I1.getMat(); Mat I1 = _I1.getMat();
...@@ -195,542 +255,539 @@ void cv::OpticalFlowDual_TVL1::operator ()(InputArray _I0, InputArray _I1, Input ...@@ -195,542 +255,539 @@ void cv::OpticalFlowDual_TVL1::operator ()(InputArray _I0, InputArray _I1, Input
merge(uxy, 2, _flow); merge(uxy, 2, _flow);
} }
namespace ////////////////////////////////////////////////////////////
{ // buildFlowMap
////////////////////////////////////////////////////////////
// buildFlowMap
struct BuildFlowMapBody : ParallelLoopBody struct BuildFlowMapBody : ParallelLoopBody
{ {
void operator() (const Range& range) const; void operator() (const Range& range) const;
Mat_<float> u1; Mat_<float> u1;
Mat_<float> u2; Mat_<float> u2;
mutable Mat_<float> map1; mutable Mat_<float> map1;
mutable Mat_<float> map2; mutable Mat_<float> 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* map1Row = map1[y];
float* map2Row = map2[y]; float* map2Row = map2[y];
for (int x = 0; x < u1.cols; ++x) for (int x = 0; x < u1.cols; ++x)
{ {
map1Row[x] = x + u1Row[x]; map1Row[x] = x + u1Row[x];
map2Row[x] = y + u2Row[x]; map2Row[x] = y + u2Row[x];
}
} }
} }
}
void buildFlowMap(const Mat_<float>& u1, const Mat_<float>& u2, Mat_<float>& map1, Mat_<float>& map2) void buildFlowMap(const Mat_<float>& u1, const Mat_<float>& u2, Mat_<float>& map1, Mat_<float>& map2)
{ {
CV_DbgAssert( u2.size() == u1.size() ); CV_DbgAssert( u2.size() == u1.size() );
CV_DbgAssert( map1.size() == u1.size() ); CV_DbgAssert( map1.size() == u1.size() );
CV_DbgAssert( map2.size() == u1.size() ); CV_DbgAssert( map2.size() == u1.size() );
BuildFlowMapBody body; BuildFlowMapBody body;
body.u1 = u1; body.u1 = u1;
body.u2 = u2; body.u2 = u2;
body.map1 = map1; body.map1 = map1;
body.map2 = map2; body.map2 = map2;
parallel_for_(Range(0, u1.rows), body); parallel_for_(Range(0, u1.rows), body);
} }
//////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////
// centeredGradient // centeredGradient
struct CenteredGradientBody : ParallelLoopBody struct CenteredGradientBody : ParallelLoopBody
{ {
void operator() (const Range& range) const; void operator() (const Range& range) const;
Mat_<float> src; Mat_<float> src;
mutable Mat_<float> dx; mutable Mat_<float> dx;
mutable Mat_<float> dy; mutable Mat_<float> dy;
}; };
void CenteredGradientBody::operator() (const Range& range) const void CenteredGradientBody::operator() (const Range& range) const
{ {
const int last_col = src.cols - 1; const int last_col = src.cols - 1;
for (int y = range.start; y < range.end; ++y) for (int y = range.start; y < range.end; ++y)
{ {
const float* srcPrevRow = src[y - 1]; const float* srcPrevRow = src[y - 1];
const float* srcCurRow = src[y]; const float* srcCurRow = src[y];
const float* srcNextRow = src[y + 1]; const float* srcNextRow = src[y + 1];
float* dxRow = dx[y]; float* dxRow = dx[y];
float* dyRow = dy[y]; float* dyRow = dy[y];
for (int x = 1; x < last_col; ++x) for (int x = 1; x < last_col; ++x)
{ {
dxRow[x] = 0.5f * (srcCurRow[x + 1] - srcCurRow[x - 1]); dxRow[x] = 0.5f * (srcCurRow[x + 1] - srcCurRow[x - 1]);
dyRow[x] = 0.5f * (srcNextRow[x] - srcPrevRow[x]); dyRow[x] = 0.5f * (srcNextRow[x] - srcPrevRow[x]);
}
} }
} }
}
void centeredGradient(const Mat_<float>& src, Mat_<float>& dx, Mat_<float>& dy) void centeredGradient(const Mat_<float>& src, Mat_<float>& dx, Mat_<float>& dy)
{ {
CV_DbgAssert( src.rows > 2 && src.cols > 2 ); CV_DbgAssert( src.rows > 2 && src.cols > 2 );
CV_DbgAssert( dx.size() == src.size() ); CV_DbgAssert( dx.size() == src.size() );
CV_DbgAssert( dy.size() == src.size() ); CV_DbgAssert( dy.size() == src.size() );
const int last_row = src.rows - 1; const int last_row = src.rows - 1;
const int last_col = src.cols - 1; const int last_col = src.cols - 1;
// compute the gradient on the center body of the image // compute the gradient on the center body of the image
{ {
CenteredGradientBody body; CenteredGradientBody body;
body.src = src; body.src = src;
body.dx = dx; body.dx = dx;
body.dy = dy; 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 // compute the gradient on the first and last rows
for (int x = 1; x < last_col; ++x) for (int x = 1; x < last_col; ++x)
{ {
dx(0, x) = 0.5f * (src(0, x + 1) - src(0, x - 1)); dx(0, x) = 0.5f * (src(0, x + 1) - src(0, x - 1));
dy(0, x) = 0.5f * (src(1, x) - src(0, x)); 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)); 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)); dy(last_row, x) = 0.5f * (src(last_row, x) - src(last_row - 1, x));
} }
// compute the gradient on the first and last columns // compute the gradient on the first and last columns
for (int y = 1; y < last_row; ++y) for (int y = 1; y < last_row; ++y)
{ {
dx(y, 0) = 0.5f * (src(y, 1) - src(y, 0)); dx(y, 0) = 0.5f * (src(y, 1) - src(y, 0));
dy(y, 0) = 0.5f * (src(y + 1, 0) - src(y - 1, 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)); 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)); dy(y, last_col) = 0.5f * (src(y + 1, last_col) - src(y - 1, last_col));
} }
// compute the gradient at the four corners // compute the gradient at the four corners
dx(0, 0) = 0.5f * (src(0, 1) - src(0, 0)); dx(0, 0) = 0.5f * (src(0, 1) - src(0, 0));
dy(0, 0) = 0.5f * (src(1, 0) - 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)); 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)); 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)); 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)); 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)); 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)); dy(last_row, last_col) = 0.5f * (src(last_row, last_col) - src(last_row - 1, last_col));
} }
//////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////
// forwardGradient // forwardGradient
struct ForwardGradientBody : ParallelLoopBody struct ForwardGradientBody : ParallelLoopBody
{ {
void operator() (const Range& range) const; void operator() (const Range& range) const;
Mat_<float> src; Mat_<float> src;
mutable Mat_<float> dx; mutable Mat_<float> dx;
mutable Mat_<float> dy; mutable Mat_<float> dy;
}; };
void ForwardGradientBody::operator() (const Range& range) const void ForwardGradientBody::operator() (const Range& range) const
{ {
const int last_col = src.cols - 1; const int last_col = src.cols - 1;
for (int y = range.start; y < range.end; ++y) for (int y = range.start; y < range.end; ++y)
{ {
const float* srcCurRow = src[y]; const float* srcCurRow = src[y];
const float* srcNextRow = src[y + 1]; const float* srcNextRow = src[y + 1];
float* dxRow = dx[y]; float* dxRow = dx[y];
float* dyRow = dy[y]; float* dyRow = dy[y];
for (int x = 0; x < last_col; ++x) for (int x = 0; x < last_col; ++x)
{ {
dxRow[x] = srcCurRow[x + 1] - srcCurRow[x]; dxRow[x] = srcCurRow[x + 1] - srcCurRow[x];
dyRow[x] = srcNextRow[x] - srcCurRow[x]; dyRow[x] = srcNextRow[x] - srcCurRow[x];
}
} }
} }
}
void forwardGradient(const Mat_<float>& src, Mat_<float>& dx, Mat_<float>& dy) void forwardGradient(const Mat_<float>& src, Mat_<float>& dx, Mat_<float>& dy)
{ {
CV_DbgAssert( src.rows > 2 && src.cols > 2 ); CV_DbgAssert( src.rows > 2 && src.cols > 2 );
CV_DbgAssert( dx.size() == src.size() ); CV_DbgAssert( dx.size() == src.size() );
CV_DbgAssert( dy.size() == src.size() ); CV_DbgAssert( dy.size() == src.size() );
const int last_row = src.rows - 1;
const int last_col = src.cols - 1;
// compute the gradient on the central body of the image const int last_row = src.rows - 1;
{ const int last_col = src.cols - 1;
ForwardGradientBody body;
body.src = src; // compute the gradient on the central body of the image
body.dx = dx; {
body.dy = dy; 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 parallel_for_(Range(0, last_row), body);
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;
}
// compute the gradient on the last column // compute the gradient on the last row
for (int y = 0; y < last_row; ++y) for (int x = 0; x < last_col; ++x)
{ {
dx(y, last_col) = 0.0f; dx(last_row, x) = src(last_row, x + 1) - src(last_row, x);
dy(y, last_col) = src(y + 1, last_col) - src(y, last_col); dy(last_row, x) = 0.0f;
} }
dx(last_row, last_col) = 0.0f; // compute the gradient on the last column
dy(last_row, last_col) = 0.0f; 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);
} }
//////////////////////////////////////////////////////////// dx(last_row, last_col) = 0.0f;
// divergence dy(last_row, last_col) = 0.0f;
}
////////////////////////////////////////////////////////////
// divergence
struct DivergenceBody : ParallelLoopBody struct DivergenceBody : ParallelLoopBody
{ {
void operator() (const Range& range) const; void operator() (const Range& range) const;
Mat_<float> v1; Mat_<float> v1;
Mat_<float> v2; Mat_<float> v2;
mutable Mat_<float> div; mutable Mat_<float> 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* v1Row = v1[y]; const float* v2CurRow = v2[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) for(int x = 1; x < v1.cols; ++x)
{ {
const float v1x = v1Row[x] - v1Row[x - 1]; const float v1x = v1Row[x] - v1Row[x - 1];
const float v2y = v2CurRow[x] - v2PrevRow[x]; const float v2y = v2CurRow[x] - v2PrevRow[x];
divRow[x] = v1x + v2y; divRow[x] = v1x + v2y;
}
} }
} }
}
void divergence(const Mat_<float>& v1, const Mat_<float>& v2, Mat_<float>& 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_<float>& v1, const Mat_<float>& v2, Mat_<float>& div)
{ {
CV_DbgAssert( v1.rows > 2 && v1.cols > 2 ); DivergenceBody body;
CV_DbgAssert( v2.size() == v1.size() );
CV_DbgAssert( div.size() == v1.size() );
{ body.v1 = v1;
DivergenceBody body; body.v2 = v2;
body.div = div;
body.v1 = v1; parallel_for_(Range(1, v1.rows), body);
body.v2 = v2; }
body.div = div;
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 // compute the divergence on the first column
for(int x = 1; x < v1.cols; ++x) for (int y = 1; y < v1.rows; ++y)
div(0, x) = v1(0, x) - v1(0, x - 1) + v2(0, x); div(y, 0) = v1(y, 0) + v2(y, 0) - v2(y - 1, 0);
// compute the divergence on the first column div(0, 0) = v1(0, 0) + v2(0, 0);
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);
}
//////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////
// calcGradRho // calcGradRho
struct CalcGradRhoBody : ParallelLoopBody struct CalcGradRhoBody : ParallelLoopBody
{ {
void operator() (const Range& range) const; void operator() (const Range& range) const;
Mat_<float> I0; Mat_<float> I0;
Mat_<float> I1w; Mat_<float> I1w;
Mat_<float> I1wx; Mat_<float> I1wx;
Mat_<float> I1wy; Mat_<float> I1wy;
Mat_<float> u1; Mat_<float> u1;
Mat_<float> u2; Mat_<float> u2;
mutable Mat_<float> grad; mutable Mat_<float> grad;
mutable Mat_<float> rho_c; mutable Mat_<float> rho_c;
}; };
void CalcGradRhoBody::operator() (const Range& range) const 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* I0Row = I0[y]; const float* I1wxRow = I1wx[y];
const float* I1wRow = I1w[y]; const float* I1wyRow = I1wy[y];
const float* I1wxRow = I1wx[y]; const float* u1Row = u1[y];
const float* I1wyRow = I1wy[y]; const float* u2Row = u2[y];
const float* u1Row = u1[y];
const float* u2Row = u2[y];
float* gradRow = grad[y]; float* gradRow = grad[y];
float* rhoRow = rho_c[y]; float* rhoRow = rho_c[y];
for (int x = 0; x < I0.cols; ++x) for (int x = 0; x < I0.cols; ++x)
{ {
const float Ix2 = I1wxRow[x] * I1wxRow[x]; const float Ix2 = I1wxRow[x] * I1wxRow[x];
const float Iy2 = I1wyRow[x] * I1wyRow[x]; const float Iy2 = I1wyRow[x] * I1wyRow[x];
// store the |Grad(I1)|^2 // store the |Grad(I1)|^2
gradRow[x] = Ix2 + Iy2; gradRow[x] = Ix2 + Iy2;
// compute the constant part of the rho function // compute the constant part of the rho function
rhoRow[x] = (I1wRow[x] - I1wxRow[x] * u1Row[x] - I1wyRow[x] * u2Row[x] - I0Row[x]); rhoRow[x] = (I1wRow[x] - I1wxRow[x] * u1Row[x] - I1wyRow[x] * u2Row[x] - I0Row[x]);
}
} }
} }
}
void calcGradRho(const Mat_<float>& I0, const Mat_<float>& I1w, const Mat_<float>& I1wx, const Mat_<float>& I1wy, const Mat_<float>& u1, const Mat_<float>& u2, void calcGradRho(const Mat_<float>& I0, const Mat_<float>& I1w, const Mat_<float>& I1wx, const Mat_<float>& I1wy, const Mat_<float>& u1, const Mat_<float>& u2,
Mat_<float>& grad, Mat_<float>& rho_c) Mat_<float>& grad, Mat_<float>& rho_c)
{ {
CV_DbgAssert( I1w.size() == I0.size() ); CV_DbgAssert( I1w.size() == I0.size() );
CV_DbgAssert( I1wx.size() == I0.size() ); CV_DbgAssert( I1wx.size() == I0.size() );
CV_DbgAssert( I1wy.size() == I0.size() ); CV_DbgAssert( I1wy.size() == I0.size() );
CV_DbgAssert( u1.size() == I0.size() ); CV_DbgAssert( u1.size() == I0.size() );
CV_DbgAssert( u2.size() == I0.size() ); CV_DbgAssert( u2.size() == I0.size() );
CV_DbgAssert( grad.size() == I0.size() ); CV_DbgAssert( grad.size() == I0.size() );
CV_DbgAssert( rho_c.size() == I0.size() ); CV_DbgAssert( rho_c.size() == I0.size() );
CalcGradRhoBody body; CalcGradRhoBody body;
body.I0 = I0; body.I0 = I0;
body.I1w = I1w; body.I1w = I1w;
body.I1wx = I1wx; body.I1wx = I1wx;
body.I1wy = I1wy; body.I1wy = I1wy;
body.u1 = u1; body.u1 = u1;
body.u2 = u2; body.u2 = u2;
body.grad = grad; body.grad = grad;
body.rho_c = rho_c; body.rho_c = rho_c;
parallel_for_(Range(0, I0.rows), body); parallel_for_(Range(0, I0.rows), body);
} }
//////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////
// estimateV // estimateV
struct EstimateVBody : ParallelLoopBody struct EstimateVBody : ParallelLoopBody
{ {
void operator() (const Range& range) const; void operator() (const Range& range) const;
Mat_<float> I1wx; Mat_<float> I1wx;
Mat_<float> I1wy; Mat_<float> I1wy;
Mat_<float> u1; Mat_<float> u1;
Mat_<float> u2; Mat_<float> u2;
Mat_<float> grad; Mat_<float> grad;
Mat_<float> rho_c; Mat_<float> rho_c;
mutable Mat_<float> v1; mutable Mat_<float> v1;
mutable Mat_<float> v2; mutable Mat_<float> v2;
float l_t; float l_t;
}; };
void EstimateVBody::operator() (const Range& range) const 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 rho = rhoRow[x] + (I1wxRow[x] * u1Row[x] + I1wyRow[x] * u2Row[x]);
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 d1 = 0.0f;
float* v2Row = v2[y]; 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]); d1 = l_t * I1wxRow[x];
d2 = l_t * I1wyRow[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<float>::epsilon())
{
float fi = -rho / gradRow[x];
d1 = fi * I1wxRow[x];
d2 = fi * I1wyRow[x];
}
v1Row[x] = u1Row[x] + d1;
v2Row[x] = u2Row[x] + d2;
} }
else if (rho > l_t * gradRow[x])
{
d1 = -l_t * I1wxRow[x];
d2 = -l_t * I1wyRow[x];
}
else if (gradRow[x] > numeric_limits<float>::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_<float>& I1wx, const Mat_<float>& I1wy, const Mat_<float>& u1, const Mat_<float>& u2, const Mat_<float>& grad, const Mat_<float>& rho_c, 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) Mat_<float>& v1, Mat_<float>& v2, float l_t)
{ {
CV_DbgAssert( I1wy.size() == I1wx.size() ); CV_DbgAssert( I1wy.size() == I1wx.size() );
CV_DbgAssert( u1.size() == I1wx.size() ); CV_DbgAssert( u1.size() == I1wx.size() );
CV_DbgAssert( u2.size() == I1wx.size() ); CV_DbgAssert( u2.size() == I1wx.size() );
CV_DbgAssert( grad.size() == I1wx.size() ); CV_DbgAssert( grad.size() == I1wx.size() );
CV_DbgAssert( rho_c.size() == I1wx.size() ); CV_DbgAssert( rho_c.size() == I1wx.size() );
CV_DbgAssert( v1.size() == I1wx.size() ); CV_DbgAssert( v1.size() == I1wx.size() );
CV_DbgAssert( v2.size() == I1wx.size() ); CV_DbgAssert( v2.size() == I1wx.size() );
EstimateVBody body; EstimateVBody body;
body.I1wx = I1wx; body.I1wx = I1wx;
body.I1wy = I1wy; body.I1wy = I1wy;
body.u1 = u1; body.u1 = u1;
body.u2 = u2; body.u2 = u2;
body.grad = grad; body.grad = grad;
body.rho_c = rho_c; body.rho_c = rho_c;
body.v1 = v1; body.v1 = v1;
body.v2 = v2; body.v2 = v2;
body.l_t = l_t; body.l_t = l_t;
parallel_for_(Range(0, I1wx.rows), body); parallel_for_(Range(0, I1wx.rows), body);
} }
//////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////
// estimateU // 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>& div_p1, const Mat_<float>& div_p2, Mat_<float>& u1, Mat_<float>& 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() ); const float* v1Row = v1[y];
CV_DbgAssert( div_p1.size() == v1.size() ); const float* v2Row = v2[y];
CV_DbgAssert( div_p2.size() == v1.size() ); const float* divP1Row = div_p1[y];
CV_DbgAssert( u1.size() == v1.size() ); const float* divP2Row = div_p2[y];
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];
float* u1Row = u1[y]; float* u1Row = u1[y];
float* u2Row = u2[y]; float* u2Row = u2[y];
for (int x = 0; x < v1.cols; ++x) for (int x = 0; x < v1.cols; ++x)
{ {
const float u1k = u1Row[x]; const float u1k = u1Row[x];
const float u2k = u2Row[x]; const float u2k = u2Row[x];
u1Row[x] = v1Row[x] + theta * divP1Row[x]; u1Row[x] = v1Row[x] + theta * divP1Row[x];
u2Row[x] = v2Row[x] + theta * divP2Row[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;
} }
//////////////////////////////////////////////////////////// return error;
// estimateDualVariables }
////////////////////////////////////////////////////////////
// estimateDualVariables
struct EstimateDualVariablesBody : ParallelLoopBody struct EstimateDualVariablesBody : ParallelLoopBody
{ {
void operator() (const Range& range) const; void operator() (const Range& range) const;
Mat_<float> u1x; Mat_<float> u1x;
Mat_<float> u1y; Mat_<float> u1y;
Mat_<float> u2x; Mat_<float> u2x;
Mat_<float> u2y; Mat_<float> u2y;
mutable Mat_<float> p11; mutable Mat_<float> p11;
mutable Mat_<float> p12; mutable Mat_<float> p12;
mutable Mat_<float> p21; mutable Mat_<float> p21;
mutable Mat_<float> p22; mutable Mat_<float> p22;
float taut; float taut;
}; };
void EstimateDualVariablesBody::operator() (const Range& range) const 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* u1xRow = u1x[y]; const float* u2xRow = u2x[y];
const float* u1yRow = u1y[y]; const float* u2yRow = u2y[y];
const float* u2xRow = u2x[y];
const float* u2yRow = u2y[y];
float* p11Row = p11[y]; float* p11Row = p11[y];
float* p12Row = p12[y]; float* p12Row = p12[y];
float* p21Row = p21[y]; float* p21Row = p21[y];
float* p22Row = p22[y]; float* p22Row = p22[y];
for (int x = 0; x < u1x.cols; ++x) for (int x = 0; x < u1x.cols; ++x)
{ {
const float g1 = static_cast<float>(hypot(u1xRow[x], u1yRow[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 g2 = static_cast<float>(hypot(u2xRow[x], u2yRow[x]));
const float ng1 = 1.0f + taut * g1; const float ng1 = 1.0f + taut * g1;
const float ng2 = 1.0f + taut * g2; const float ng2 = 1.0f + taut * g2;
p11Row[x] = (p11Row[x] + taut * u1xRow[x]) / ng1; p11Row[x] = (p11Row[x] + taut * u1xRow[x]) / ng1;
p12Row[x] = (p12Row[x] + taut * u1yRow[x]) / ng1; p12Row[x] = (p12Row[x] + taut * u1yRow[x]) / ng1;
p21Row[x] = (p21Row[x] + taut * u2xRow[x]) / ng2; p21Row[x] = (p21Row[x] + taut * u2xRow[x]) / ng2;
p22Row[x] = (p22Row[x] + taut * u2yRow[x]) / ng2; p22Row[x] = (p22Row[x] + taut * u2yRow[x]) / ng2;
}
} }
} }
}
void estimateDualVariables(const Mat_<float>& u1x, const Mat_<float>& u1y, const Mat_<float>& u2x, const Mat_<float>& u2y, 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) Mat_<float>& p11, Mat_<float>& p12, Mat_<float>& p21, Mat_<float>& p22, float taut)
{ {
CV_DbgAssert( u1y.size() == u1x.size() ); CV_DbgAssert( u1y.size() == u1x.size() );
CV_DbgAssert( u2x.size() == u1x.size() ); CV_DbgAssert( u2x.size() == u1x.size() );
CV_DbgAssert( u2y.size() == u1x.size() ); CV_DbgAssert( u2y.size() == u1x.size() );
CV_DbgAssert( p11.size() == u1x.size() ); CV_DbgAssert( p11.size() == u1x.size() );
CV_DbgAssert( p12.size() == u1x.size() ); CV_DbgAssert( p12.size() == u1x.size() );
CV_DbgAssert( p21.size() == u1x.size() ); CV_DbgAssert( p21.size() == u1x.size() );
CV_DbgAssert( p22.size() == u1x.size() ); CV_DbgAssert( p22.size() == u1x.size() );
EstimateDualVariablesBody body; EstimateDualVariablesBody body;
body.u1x = u1x; body.u1x = u1x;
body.u1y = u1y; body.u1y = u1y;
body.u2x = u2x; body.u2x = u2x;
body.u2y = u2y; body.u2y = u2y;
body.p11 = p11; body.p11 = p11;
body.p12 = p12; body.p12 = p12;
body.p21 = p21; body.p21 = p21;
body.p22 = p22; body.p22 = p22;
body.taut = taut; body.taut = taut;
parallel_for_(Range(0, u1x.rows), body); parallel_for_(Range(0, u1x.rows), body);
}
} }
void cv::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)
{ {
const float scaledEpsilon = static_cast<float>(epsilon * epsilon * I0.size().area()); const float scaledEpsilon = static_cast<float>(epsilon * epsilon * I0.size().area());
...@@ -818,21 +875,12 @@ void cv::OpticalFlowDual_TVL1::procOneScale(const Mat_<float>& I0, const Mat_<fl ...@@ -818,21 +875,12 @@ void cv::OpticalFlowDual_TVL1::procOneScale(const Mat_<float>& I0, const Mat_<fl
} }
} }
namespace void OpticalFlowDual_TVL1::collectGarbage()
{
template <typename T> void releaseVector(vector<T>& v)
{
vector<T> empty;
empty.swap(v);
}
}
void cv::OpticalFlowDual_TVL1::collectGarbage()
{ {
releaseVector(I0s); I0s.clear();
releaseVector(I1s); I1s.clear();
releaseVector(u1s); u1s.clear();
releaseVector(u2s); u2s.clear();
I1x_buf.release(); I1x_buf.release();
I1y_buf.release(); I1y_buf.release();
...@@ -863,3 +911,27 @@ void cv::OpticalFlowDual_TVL1::collectGarbage() ...@@ -863,3 +911,27 @@ void cv::OpticalFlowDual_TVL1::collectGarbage()
u2x_buf.release(); u2x_buf.release();
u2y_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<DenseOpticalFlow> cv::createOptFlow_DualTVL1()
{
return new OpticalFlowDual_TVL1;
}
...@@ -152,9 +152,9 @@ TEST(Video_calcOpticalFlowDual_TVL1, Regression) ...@@ -152,9 +152,9 @@ TEST(Video_calcOpticalFlowDual_TVL1, Regression)
ASSERT_FALSE(frame2.empty()); ASSERT_FALSE(frame2.empty());
Mat_<Point2f> flow; Mat_<Point2f> flow;
OpticalFlowDual_TVL1 tvl1; Ptr<DenseOpticalFlow> tvl1 = createOptFlow_DualTVL1();
tvl1(frame1, frame2, flow); tvl1->calc(frame1, frame2, flow);
#ifdef DUMP #ifdef DUMP
writeOpticalFlowToFile(flow, gold_flow_path); writeOpticalFlowToFile(flow, gold_flow_path);
......
...@@ -173,10 +173,10 @@ int main(int argc, const char* argv[]) ...@@ -173,10 +173,10 @@ int main(int argc, const char* argv[])
} }
Mat_<Point2f> flow; Mat_<Point2f> flow;
OpticalFlowDual_TVL1 tvl1; Ptr<DenseOpticalFlow> tvl1 = createOptFlow_DualTVL1();
const double start = (double)getTickCount(); const double start = (double)getTickCount();
tvl1(frame0, frame1, flow); tvl1->calc(frame0, frame1, flow);
const double timeSec = (getTickCount() - start) / getTickFrequency(); const double timeSec = (getTickCount() - start) / getTickFrequency();
cout << "calcOpticalFlowDual_TVL1 : " << timeSec << " sec" << endl; cout << "calcOpticalFlowDual_TVL1 : " << timeSec << " sec" << endl;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment