Skip to content

O
opencv
Commit 44e04da7 authored by Vadim Pisarevsky's avatar Vadim Pisarevsky
Browse files
Options

Merge pull request #2942 from ernest-galbrun:tvl1_chambolle

parents b2cd954f 2f077fcd
Hide whitespace changes
Inline Side-by-side
Showing with 216 additions and 59 deletions
modules/cudaoptflow/include/opencv2/cudaoptflow.hpp
View file @ 44e04da7
......@@ -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;
......
modules/cudaoptflow/perf/perf_optflow.cpp
View file @ 44e04da7
......@@ -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);
......
modules/cudaoptflow/src/cuda/tvl1flow.cu
View file @ 44e04da7
......@@ -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() );
......
modules/cudaoptflow/src/tvl1flow.cpp
View file @ 44e04da7
......@@ -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();
}
......
modules/cudaoptflow/test/test_optflow.cpp
View file @ 44e04da7
......@@ -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);
......
modules/cudastereo/src/cuda/stereocsbp.cu
View file @ 44e04da7
......@@ -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>
......
modules/video/src/tvl1flow.cpp
View file @ 44e04da7
......@@ -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
......
modules/video/test/test_tvl1optflow.cpp
View file @ 44e04da7
......@@ -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
}
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