__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,
- const PtrStepf p31, const PtrStepf p32,
- PtrStepf u1, PtrStepf u2, PtrStepf u3, PtrStepf error,
+ 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 float I1wyVal = I1wy(y, x);
const float gradVal = grad(y, x);
const float u1OldVal = u1(y, x);
- const float u2OldVal = u2(y, x);
- const float u3OldVal = u3(y, x);
+ const float u2OldVal = u2(y, x);
+ const float u3OldVal = u3(y, x);
const float rho = rho_c(y, x) + (I1wxVal * u1OldVal + I1wyVal * u2OldVal + gamma * u3OldVal);
float d1 = 0.0f;
float d2 = 0.0f;
- float d3 = 0.0f;
+ float d3 = 0.0f;
if (rho < -l_t * gradVal)
{
d1 = l_t * I1wxVal;
d2 = l_t * I1wyVal;
- d3 = l_t * gamma;
+ d3 = l_t * gamma;
}
else if (rho > l_t * gradVal)
{
d1 = -l_t * I1wxVal;
d2 = -l_t * I1wyVal;
- d3 = -l_t * gamma;
+ d3 = -l_t * gamma;
}
else if (gradVal > numeric_limits<float>::epsilon())
{
const float fi = -rho / gradVal;
d1 = fi * I1wxVal;
d2 = fi * I1wyVal;
- d3 = fi * gamma;
+ d3 = fi * gamma;
}
const float v1 = u1OldVal + d1;
- const float v2 = u2OldVal + d2;
- const float v3 = u3OldVal + d3;
+ 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 = divergence(p31, p32, y, x);
+ const float div_p2 = divergence(p21, p22, y, x);
+ const float div_p3 = divergence(p31, p32, y, x);
// 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 = v3 + theta * div_p3;
+ const float u2NewVal = v2 + theta * div_p2;
+ const float u3NewVal = v3 + theta * div_p3;
u1(y, x) = u1NewVal;
- u2(y, x) = u2NewVal;
- u3(y, x) = u3NewVal;
+ u2(y, x) = u2NewVal;
+ u3(y, x) = u3NewVal;
if (calcError)
{
const float n1 = (u1OldVal - u1NewVal) * (u1OldVal - u1NewVal);
- const float n2 = (u2OldVal - u2NewVal) * (u2OldVal - u2NewVal);
- const float n3 = 0;// (u3OldVal - u3NewVal) * (u3OldVal - u3NewVal);
+ const float n2 = (u2OldVal - u2NewVal) * (u2OldVal - u2NewVal);
+ const float n3 = 0;// (u3OldVal - u3NewVal) * (u3OldVal - u3NewVal);
error(y, x) = n1 + n2 + n3;
}
}
void estimateU(PtrStepSzf I1wx, PtrStepSzf I1wy,
PtrStepSzf grad, PtrStepSzf rho_c,
- PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32,
- PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf error,
+ 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);
namespace tvl1flow
{
- __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)
+ __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 int x = blockIdx.x * blockDim.x + threadIdx.x;
const int y = blockIdx.y * blockDim.y + threadIdx.y;
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 = u3(y, ::min(x + 1, u1.cols - 1)) - u3(y, x);
- const float u3y = u3(::min(y + 1, u1.rows - 1), x) - u3(y, x);
+ const float u3x = u3(y, ::min(x + 1, u1.cols - 1)) - u3(y, x);
+ const float u3y = u3(::min(y + 1, u1.rows - 1), x) - u3(y, x);
const float g1 = ::hypotf(u1x, u1y);
- const float g2 = ::hypotf(u2x, u2y);
- const float g3 = ::hypotf(u3x, u3y);
+ const float g2 = ::hypotf(u2x, u2y);
+ const float g3 = ::hypotf(u3x, u3y);
const float ng1 = 1.0f + taut * g1;
- const float ng2 = 1.0f + taut * g2;
- const float ng3 = 1.0f + taut * g3;
+ const float ng2 = 1.0f + taut * g2;
+ const float ng3 = 1.0f + taut * g3;
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;
- p31(y, x) = (p31(y, x) + taut * u3x) / ng3;
- p32(y, x) = (p32(y, x) + taut * u3y) / ng3;
+ p22(y, x) = (p22(y, x) + taut * u2y) / ng2;
+ 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 u3, PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32, float taut)
+ void estimateDualVariables(PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32, float taut)
{
const dim3 block(32, 8);
const dim3 grid(divUp(u1.cols, block.x), divUp(u1.rows, block.y));
{
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)
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 = use_gamma ? rhoRow[x] + (I1wxRow[x] * u1Row[x] + I1wyRow[x] * u2Row[x]) + gamma * u3Row[x] :
////////////////////////////////////////////////////////////
// estimateU
-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,
+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)
{
u1Row[x] = v1Row[x] + theta * divP1Row[x];
u2Row[x] = v2Row[x] + theta * divP2Row[x];
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);
}
}
}
-void estimateDualVariables(const Mat_<float>& u1x, const Mat_<float>& u1y,
+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>& p11, Mat_<float>& p12,
+ Mat_<float>& p21, Mat_<float>& p22,
Mat_<float>& p31, Mat_<float>& p32,
float taut, bool use_gamma)
{