added patch error calculation to gpu::PyrLKOpticalFlow

diff --git a/modules/gpu/include/opencv2/gpu/gpu.hpp b/modules/gpu/include/opencv2/gpu/gpu.hpp
index b389d6f..cc6d607 100644 (file)
--- a/modules/gpu/include/opencv2/gpu/gpu.hpp
+++ b/modules/gpu/include/opencv2/gpu/gpu.hpp
@@ -1817,6 +1817,7 @@ public:
         derivLambda = 0.5;\r
         useInitialFlow = false;\r
         minEigThreshold = 1e-4f;\r
+        getMinEigenVals = false;\r
     }\r
 \r
     void sparse(const GpuMat& prevImg, const GpuMat& nextImg, const GpuMat& prevPts, GpuMat& nextPts,\r
@@ -1830,6 +1831,7 @@ public:
     double derivLambda;\r
     bool useInitialFlow;\r
     float minEigThreshold;\r
+    bool getMinEigenVals;\r
 \r
     void releaseMemory()\r
     {\r

diff --git a/modules/gpu/src/cuda/pyrlk.cu b/modules/gpu/src/cuda/pyrlk.cu
index 32dabb1..e41e4ca 100644 (file)
--- a/modules/gpu/src/cuda/pyrlk.cu
+++ b/modules/gpu/src/cuda/pyrlk.cu
@@ -274,9 +274,39 @@ namespace cv { namespace gpu { namespace device
             }\r
         }\r
 \r
+        __device__ void reduce(float& val1, float* smem1, int tid)\r
+        {\r
+            smem1[tid] = val1;\r
+            __syncthreads();\r
+\r
+            if (tid < 128) \r
+            { \r
+                smem1[tid] = val1 += smem1[tid + 128]; \r
+            } \r
+            __syncthreads();\r
+\r
+            if (tid < 64) \r
+            { \r
+                smem1[tid] = val1 += smem1[tid + 64]; \r
+            } \r
+            __syncthreads();\r
+\r
+            if (tid < 32)\r
+            {\r
+                volatile float* vmem1 = smem1;\r
+\r
+                vmem1[tid] = val1 += vmem1[tid + 32]; \r
+                vmem1[tid] = val1 += vmem1[tid + 16]; \r
+                vmem1[tid] = val1 += vmem1[tid + 8]; \r
+                vmem1[tid] = val1 += vmem1[tid + 4];\r
+                vmem1[tid] = val1 += vmem1[tid + 2]; \r
+                vmem1[tid] = val1 += vmem1[tid + 1]; \r
+            }\r
+        }\r
+\r
         #define SCALE (1.0f / (1 << 20))\r
 \r
-        template <int PATCH_X, int PATCH_Y, bool calcErr>\r
+        template <int PATCH_X, int PATCH_Y, bool calcErr, bool GET_MIN_EIGENVALS>\r
         __global__ void lkSparse(const PtrStepb I, const PtrStepb J, const PtrStep<short> dIdx, const PtrStep<short> dIdy,\r
             const float2* prevPts, float2* nextPts, uchar* status, float* err, const int level, const int rows, const int cols)\r
         {\r
@@ -349,7 +379,7 @@ namespace cv { namespace gpu { namespace device
                 float D = A11 * A22 - A12 * A12;\r
                 float minEig = (A22 + A11 - ::sqrtf((A11 - A22) * (A11 - A22) + 4.f * A12 * A12)) / (2 * c_winSize_x * c_winSize_y);\r
             \r
-                if (calcErr && tid == 0) \r
+                if (calcErr && GET_MIN_EIGENVALS && tid == 0) \r
                     err[blockIdx.x] = minEig;\r
 \r
                 if (minEig < c_minEigThreshold || D < numeric_limits<float>::epsilon())\r
@@ -377,7 +407,7 @@ namespace cv { namespace gpu { namespace device
             bool status_ = true;\r
 \r
             for (int k = 0; k < c_iters; ++k)\r
-            {                \r
+            {\r
                 if (nextPt.x < -c_winSize_x || nextPt.x >= cols || nextPt.y < -c_winSize_y || nextPt.y >= rows)\r
                 {\r
                     status_ = false;\r
@@ -415,38 +445,76 @@ namespace cv { namespace gpu { namespace device
                 nextPt.y += delta.y;\r
 \r
                 if (::fabs(delta.x) < 0.01f && ::fabs(delta.y) < 0.01f)\r
+                {\r
+                    nextPt.x -= delta.x * 0.5f;\r
+                    nextPt.y -= delta.y * 0.5f;\r
                     break;\r
+                }\r
             }\r
 \r
-            if (tid == 0)\r
+            if (nextPt.x < -c_winSize_x || nextPt.x >= cols || nextPt.y < -c_winSize_y || nextPt.y >= rows)\r
+                status_ = false;\r
+\r
+            // TODO : Why do we compute patch error in shifted window?\r
+            nextPt.x += c_halfWin_x;\r
+            nextPt.y += c_halfWin_y;\r
+\r
+            float errval = 0.f;\r
+            if (calcErr && !GET_MIN_EIGENVALS && status_)\r
             {\r
-                nextPt.x += c_halfWin_x;\r
-                nextPt.y += c_halfWin_y;\r
+                for (int y = threadIdx.y, i = 0; y < c_winSize_y; y += blockDim.y, ++i)\r
+                {\r
+                    for (int x = threadIdx.x, j = 0; x < c_winSize_x_cn; x += blockDim.x, ++j)\r
+                    {\r
+                        int diff = linearFilter(J, nextPt, x, y) - I_patch[i][j];\r
+                        errval += ::fabsf((float)diff);\r
+                    }\r
+                }\r
 \r
-                nextPts[blockIdx.x] = nextPt;\r
+                reduce(errval, smem1, tid);\r
+\r
+                errval /= 32 * c_winSize_x_cn * c_winSize_y;\r
+            }\r
+\r
+            if (tid == 0)\r
+            {\r
                 status[blockIdx.x] = status_;\r
+                nextPts[blockIdx.x] = nextPt;\r
+\r
+                if (calcErr && !GET_MIN_EIGENVALS)\r
+                    err[blockIdx.x] = errval;\r
             }\r
         }\r
 \r
         template <int PATCH_X, int PATCH_Y>\r
         void lkSparse_caller(DevMem2Db I, DevMem2Db J, DevMem2D_<short> dIdx, DevMem2D_<short> dIdy,\r
-            const float2* prevPts, float2* nextPts, uchar* status, float* err, int ptcount, \r
+            const float2* prevPts, float2* nextPts, uchar* status, float* err, bool GET_MIN_EIGENVALS, int ptcount, \r
             int level, dim3 block, cudaStream_t stream)\r
         {\r
             dim3 grid(ptcount);\r
 \r
-            if (err)\r
+            if (level == 0 && err)\r
             {\r
-                cudaSafeCall( cudaFuncSetCacheConfig(lkSparse<PATCH_X, PATCH_Y, true>, cudaFuncCachePreferL1) );\r
+                if (GET_MIN_EIGENVALS)\r
+                {\r
+                    cudaSafeCall( cudaFuncSetCacheConfig(lkSparse<PATCH_X, PATCH_Y, true, true>, cudaFuncCachePreferL1) );\r
+\r
+                    lkSparse<PATCH_X, PATCH_Y, true, true><<<grid, block>>>(I, J, dIdx, dIdy,\r
+                        prevPts, nextPts, status, err, level, I.rows, I.cols);\r
+                }\r
+                else\r
+                {\r
+                    cudaSafeCall( cudaFuncSetCacheConfig(lkSparse<PATCH_X, PATCH_Y, true, false>, cudaFuncCachePreferL1) );\r
 \r
-                lkSparse<PATCH_X, PATCH_Y, true><<<grid, block>>>(I, J, dIdx, dIdy,\r
-                    prevPts, nextPts, status, err, level, I.rows, I.cols);\r
+                    lkSparse<PATCH_X, PATCH_Y, true, false><<<grid, block>>>(I, J, dIdx, dIdy,\r
+                        prevPts, nextPts, status, err, level, I.rows, I.cols);\r
+                }\r
             }\r
             else\r
             {\r
-                cudaSafeCall( cudaFuncSetCacheConfig(lkSparse<PATCH_X, PATCH_Y, false>, cudaFuncCachePreferL1) );\r
+                cudaSafeCall( cudaFuncSetCacheConfig(lkSparse<PATCH_X, PATCH_Y, false, false>, cudaFuncCachePreferL1) );\r
 \r
-                lkSparse<PATCH_X, PATCH_Y, false><<<grid, block>>>(I, J, dIdx, dIdy,\r
+                lkSparse<PATCH_X, PATCH_Y, false, false><<<grid, block>>>(I, J, dIdx, dIdy,\r
                         prevPts, nextPts, status, err, level, I.rows, I.cols);\r
             }\r
 \r
@@ -457,11 +525,11 @@ namespace cv { namespace gpu { namespace device
         }\r
 \r
         void lkSparse_gpu(DevMem2Db I, DevMem2Db J, DevMem2D_<short> dIdx, DevMem2D_<short> dIdy,\r
-            const float2* prevPts, float2* nextPts, uchar* status, float* err, int ptcount, \r
+            const float2* prevPts, float2* nextPts, uchar* status, float* err, bool GET_MIN_EIGENVALS, int ptcount, \r
             int level, dim3 block, dim3 patch, cudaStream_t stream)\r
         {\r
             typedef void (*func_t)(DevMem2Db I, DevMem2Db J, DevMem2D_<short> dIdx, DevMem2D_<short> dIdy,\r
-                const float2* prevPts, float2* nextPts, uchar* status, float* err, int ptcount, \r
+                const float2* prevPts, float2* nextPts, uchar* status, float* err, bool GET_MIN_EIGENVALS, int ptcount, \r
                 int level, dim3 block, cudaStream_t stream);\r
 \r
             static const func_t funcs[5][5] = \r
@@ -474,11 +542,11 @@ namespace cv { namespace gpu { namespace device
             };            \r
 \r
             funcs[patch.y - 1][patch.x - 1](I, J, dIdx, dIdy,\r
-                prevPts, nextPts, status, err, ptcount, \r
+                prevPts, nextPts, status, err, GET_MIN_EIGENVALS, ptcount, \r
                 level, block, stream);\r
         }\r
 \r
-        template <bool calcErr>\r
+        template <bool calcErr, bool GET_MIN_EIGENVALS>\r
         __global__ void lkDense(const PtrStepb I, const PtrStepb J, const PtrStep<short> dIdx, const PtrStep<short> dIdy,\r
             PtrStepf u, PtrStepf v, PtrStepf err, const int rows, const int cols)\r
         {\r
@@ -515,7 +583,7 @@ namespace cv { namespace gpu { namespace device
                 float D = A11 * A22 - A12 * A12;\r
                 float minEig = (A22 + A11 - ::sqrtf((A11 - A22) * (A11 - A22) + 4.f * A12 * A12)) / (2 * c_winSize_x * c_winSize_y);\r
 \r
-                if (calcErr)\r
+                if (calcErr && GET_MIN_EIGENVALS)\r
                     err(y, x) = minEig;\r
             \r
                 if (minEig < c_minEigThreshold || D < numeric_limits<float>::epsilon())\r
@@ -565,30 +633,63 @@ namespace cv { namespace gpu { namespace device
 \r
                 if (::fabs(delta.x) < 0.01f && ::fabs(delta.y) < 0.01f)\r
                     break;\r
-            }\r
+            }            \r
+\r
+            // TODO : Why do we compute patch error in shifted window?\r
+            nextPt.x += c_halfWin_x;\r
+            nextPt.y += c_halfWin_y;\r
 \r
-            u(y, x) = nextPt.x - x + c_halfWin_x;\r
-            v(y, x) = nextPt.y - y + c_halfWin_y;\r
+            u(y, x) = nextPt.x - x;\r
+            v(y, x) = nextPt.y - y;            \r
+\r
+            if (calcErr && !GET_MIN_EIGENVALS)\r
+            {\r
+                float errval = 0.0f;\r
+\r
+                for (int i = 0; i < c_winSize_y; ++i)\r
+                {\r
+                    for (int j = 0; j < c_winSize_x; ++j)\r
+                    {\r
+                        int I_val = I(y - c_halfWin_y + i, x - c_halfWin_x + j);\r
+                        int diff = linearFilter(J, nextPt, j, i) - CV_DESCALE(I_val * (1 << W_BITS), W_BITS1 - 5);\r
+                        errval += ::fabsf((float)diff);\r
+                    }\r
+                }\r
+\r
+                errval /= 32 * c_winSize_x_cn * c_winSize_y;\r
+\r
+                err(y, x) = errval;\r
+            }\r
         }\r
 \r
         void lkDense_gpu(DevMem2Db I, DevMem2Db J, DevMem2D_<short> dIdx, DevMem2D_<short> dIdy, \r
-            DevMem2Df u, DevMem2Df v, DevMem2Df* err, cudaStream_t stream)\r
+            DevMem2Df u, DevMem2Df v, DevMem2Df* err, bool GET_MIN_EIGENVALS, cudaStream_t stream)\r
         {\r
             dim3 block(32, 8);\r
             dim3 grid(divUp(I.cols, block.x), divUp(I.rows, block.y));\r
 \r
             if (err)\r
             {\r
-                cudaSafeCall( cudaFuncSetCacheConfig(lkDense<true>, cudaFuncCachePreferL1) );\r
+                if (GET_MIN_EIGENVALS)\r
+                {\r
+                    cudaSafeCall( cudaFuncSetCacheConfig(lkDense<true, true>, cudaFuncCachePreferL1) );\r
 \r
-                lkDense<true><<<grid, block, 0, stream>>>(I, J, dIdx, dIdy, u, v, *err, I.rows, I.cols);\r
-                cudaSafeCall( cudaGetLastError() );\r
+                    lkDense<true, true><<<grid, block, 0, stream>>>(I, J, dIdx, dIdy, u, v, *err, I.rows, I.cols);\r
+                    cudaSafeCall( cudaGetLastError() );\r
+                }\r
+                else\r
+                {\r
+                    cudaSafeCall( cudaFuncSetCacheConfig(lkDense<true, false>, cudaFuncCachePreferL1) );\r
+\r
+                    lkDense<true, false><<<grid, block, 0, stream>>>(I, J, dIdx, dIdy, u, v, *err, I.rows, I.cols);\r
+                    cudaSafeCall( cudaGetLastError() );\r
+                }\r
             }\r
             else\r
             {\r
-                cudaSafeCall( cudaFuncSetCacheConfig(lkDense<false>, cudaFuncCachePreferL1) );\r
+                cudaSafeCall( cudaFuncSetCacheConfig(lkDense<false, false>, cudaFuncCachePreferL1) );\r
 \r
-                lkDense<false><<<grid, block, 0, stream>>>(I, J, dIdx, dIdy, u, v, PtrStepf(), I.rows, I.cols);\r
+                lkDense<false, false><<<grid, block, 0, stream>>>(I, J, dIdx, dIdy, u, v, PtrStepf(), I.rows, I.cols);\r
                 cudaSafeCall( cudaGetLastError() );\r
             }\r
 \r

diff --git a/modules/gpu/src/pyrlk.cpp b/modules/gpu/src/pyrlk.cpp
index 37427c5..83fc1f7 100644 (file)
--- a/modules/gpu/src/pyrlk.cpp
+++ b/modules/gpu/src/pyrlk.cpp
@@ -63,11 +63,11 @@ namespace cv { namespace gpu { namespace device
             cudaStream_t stream = 0);\r
 \r
         void lkSparse_gpu(DevMem2Db I, DevMem2Db J, DevMem2D_<short> dIdx, DevMem2D_<short> dIdy,\r
-            const float2* prevPts, float2* nextPts, uchar* status, float* err, int ptcount, \r
+            const float2* prevPts, float2* nextPts, uchar* status, float* err, bool GET_MIN_EIGENVALS, int ptcount, \r
             int level, dim3 block, dim3 patch, cudaStream_t stream = 0);\r
 \r
         void lkDense_gpu(DevMem2Db I, DevMem2Db J, DevMem2D_<short> dIdx, DevMem2D_<short> dIdy, \r
-            DevMem2Df u, DevMem2Df v, DevMem2Df* err, cudaStream_t stream = 0);\r
+            DevMem2Df u, DevMem2Df v, DevMem2Df* err, bool GET_MIN_EIGENVALS, cudaStream_t stream = 0);\r
     }\r
 }}}\r
 \r
@@ -205,7 +205,7 @@ void cv::gpu::PyrLKOpticalFlow::sparse(const GpuMat& prevImg, const GpuMat& next
         calcSharrDeriv(prevPyr_[level], dIdx, dIdy);\r
 \r
         lkSparse_gpu(prevPyr_[level], nextPyr_[level], dIdx, dIdy, \r
-            prevPts.ptr<float2>(), nextPts.ptr<float2>(), status.ptr(), level == 0 && err ? err->ptr<float>() : 0, prevPts.cols, \r
+            prevPts.ptr<float2>(), nextPts.ptr<float2>(), status.ptr(), level == 0 && err ? err->ptr<float>() : 0, getMinEigenVals, prevPts.cols, \r
             level, block, patch);\r
     }\r
 }\r
@@ -272,7 +272,7 @@ void cv::gpu::PyrLKOpticalFlow::dense(const GpuMat& prevImg, const GpuMat& nextI
         calcSharrDeriv(prevPyr_[level], dIdx, dIdy);\r
 \r
         lkDense_gpu(prevPyr_[level], nextPyr_[level], dIdx, dIdy, uPyr_[level], vPyr_[level], \r
-            level == 0 && err ? &derr : 0);\r
+            level == 0 && err ? &derr : 0, getMinEigenVals);\r
 \r
         if (level == 0)\r
         {\r

diff --git a/modules/gpu/test/test_video.cpp b/modules/gpu/test/test_video.cpp
index a63c4fd..fa98d1e 100644 (file)
--- a/modules/gpu/test/test_video.cpp
+++ b/modules/gpu/test/test_video.cpp
@@ -358,7 +358,7 @@ PARAM_TEST_CASE(PyrLKOpticalFlowSparse, cv::gpu::DeviceInfo, bool)
         cv::goodFeaturesToTrack(gray_frame, pts, 1000, 0.01, 0.0);\r
 \r
         cv::calcOpticalFlowPyrLK(frame0, frame1, pts, nextPts_gold, status_gold, err_gold, cv::Size(21, 21), 3, \r
-            cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 30, 0.01), 0.5, CV_LKFLOW_GET_MIN_EIGENVALS);\r
+            cv::TermCriteria(cv::TermCriteria::COUNT + cv::TermCriteria::EPS, 30, 0.01), 0.5);\r
     }\r
 };\r
 \r
@@ -410,7 +410,7 @@ TEST_P(PyrLKOpticalFlowSparse, Accuracy)
             bool eq = std::abs(a.x - b.x) < 1 && std::abs(a.y - b.y) < 1;\r
             float errdiff = std::abs(err[i] - err_gold[i]);\r
 \r
-            if (!eq || errdiff > 1e-4)\r
+            if (!eq || errdiff > 1e-1)\r
                 ++mistmatch;\r
         }\r
     }\r

diff --git a/samples/gpu/performance/tests.cpp b/samples/gpu/performance/tests.cpp
index cb51abd..d89b7d1 100644 (file)
--- a/samples/gpu/performance/tests.cpp
+++ b/samples/gpu/performance/tests.cpp
@@ -1157,10 +1157,12 @@ TEST(PyrLKOpticalFlow)
         vector<Point2f> nextPts;\r
         vector<unsigned char> status;\r
 \r
-        calcOpticalFlowPyrLK(frame0, frame1, pts, nextPts, status, noArray());\r
+        vector<float> err;\r
+\r
+        calcOpticalFlowPyrLK(frame0, frame1, pts, nextPts, status, err);\r
 \r
         CPU_ON;\r
-        calcOpticalFlowPyrLK(frame0, frame1, pts, nextPts, status, noArray());\r
+        calcOpticalFlowPyrLK(frame0, frame1, pts, nextPts, status, err);\r
         CPU_OFF;\r
 \r
         gpu::PyrLKOpticalFlow d_pyrLK;\r
@@ -1176,10 +1178,10 @@ TEST(PyrLKOpticalFlow)
         gpu::GpuMat d_status;\r
         gpu::GpuMat d_err;\r
 \r
-        d_pyrLK.sparse(d_frame0, d_frame1, d_pts, d_nextPts, d_status);\r
+        d_pyrLK.sparse(d_frame0, d_frame1, d_pts, d_nextPts, d_status, &d_err);\r
 \r
         GPU_ON;\r
-        d_pyrLK.sparse(d_frame0, d_frame1, d_pts, d_nextPts, d_status);\r
+        d_pyrLK.sparse(d_frame0, d_frame1, d_pts, d_nextPts, d_status, &d_err);\r
         GPU_OFF;\r
     }\r
 }\r