diff --git a/modules/ocl/include/opencv2/ocl/ocl.hpp b/modules/ocl/include/opencv2/ocl/ocl.hpp index 7395c7bcc2..aecf2be748 100644 --- a/modules/ocl/include/opencv2/ocl/ocl.hpp +++ b/modules/ocl/include/opencv2/ocl/ocl.hpp @@ -1410,6 +1410,69 @@ namespace cv oclMat vPyr_[2]; bool isDeviceArch11_; }; + + class CV_EXPORTS FarnebackOpticalFlow + { + public: + FarnebackOpticalFlow() + { + numLevels = 5; + pyrScale = 0.5; + fastPyramids = false; + winSize = 13; + numIters = 10; + polyN = 5; + polySigma = 1.1; + flags = 0; + } + + int numLevels; + double pyrScale; + bool fastPyramids; + int winSize; + int numIters; + int polyN; + double polySigma; + int flags; + + void operator ()(const oclMat &frame0, const oclMat &frame1, oclMat &flowx, oclMat &flowy); + + void releaseMemory() + { + frames_[0].release(); + frames_[1].release(); + pyrLevel_[0].release(); + pyrLevel_[1].release(); + M_.release(); + bufM_.release(); + R_[0].release(); + R_[1].release(); + blurredFrame_[0].release(); + blurredFrame_[1].release(); + pyramid0_.clear(); + pyramid1_.clear(); + } + + private: + void prepareGaussian( + int n, double sigma, float *g, float *xg, float *xxg, + double &ig11, double &ig03, double &ig33, double &ig55); + + void setPolynomialExpansionConsts(int n, double sigma); + + void updateFlow_boxFilter( + const oclMat& R0, const oclMat& R1, oclMat& flowx, oclMat &flowy, + oclMat& M, oclMat &bufM, int blockSize, bool updateMatrices); + + void updateFlow_gaussianBlur( + const oclMat& R0, const oclMat& R1, oclMat& flowx, oclMat& flowy, + oclMat& M, oclMat &bufM, int blockSize, bool updateMatrices); + + oclMat frames_[2]; + oclMat pyrLevel_[2], M_, bufM_, R_[2], blurredFrame_[2]; + std::vector pyramid0_, pyramid1_; + }; + //////////////// build warping maps //////////////////// //! builds plane warping maps CV_EXPORTS void buildWarpPlaneMaps(Size src_size, Rect dst_roi, const Mat &K, const Mat &R, const Mat &T, float scale, oclMat &map_x, oclMat &map_y); diff --git a/modules/ocl/perf/perf_opticalflow.cpp b/modules/ocl/perf/perf_opticalflow.cpp index 97283b206c..4b987f4817 100644 --- a/modules/ocl/perf/perf_opticalflow.cpp +++ b/modules/ocl/perf/perf_opticalflow.cpp @@ -225,4 +225,133 @@ PERFTEST(tvl1flow) TestSystem::instance().ExceptedMatSimilar(gold[0], flowx, 3e-3); TestSystem::instance().ExceptedMatSimilar(gold[1], flowy, 3e-3); -} \ No newline at end of file +} + +///////////// FarnebackOpticalFlow //////////////////////// +PERFTEST(FarnebackOpticalFlow) +{ + cv::Mat frame0 = imread("rubberwhale1.png", cv::IMREAD_GRAYSCALE); + ASSERT_FALSE(frame0.empty()); + + cv::Mat frame1 = imread("rubberwhale2.png", cv::IMREAD_GRAYSCALE); + ASSERT_FALSE(frame1.empty()); + + cv::ocl::oclMat d_frame0(frame0), d_frame1(frame1); + + int polyNs[2] = { 5, 7 }; + double polySigmas[2] = { 1.1, 1.5 }; + int farneFlags[2] = { 0, cv::OPTFLOW_FARNEBACK_GAUSSIAN }; + bool UseInitFlows[2] = { false, true }; + double pyrScale = 0.5; + + string farneFlagStrs[2] = { "BoxFilter", "GaussianBlur" }; + string useInitFlowStrs[2] = { "", "UseInitFlow" }; + + for ( int i = 0; i < 2; ++i) + { + int polyN = polyNs[i]; + double polySigma = polySigmas[i]; + + for ( int j = 0; j < 2; ++j) + { + int flags = farneFlags[j]; + + for ( int k = 0; k < 2; ++k) + { + bool useInitFlow = UseInitFlows[k]; + SUBTEST << "polyN(" << polyN << "); " << farneFlagStrs[j] << "; " << useInitFlowStrs[k]; + + cv::ocl::FarnebackOpticalFlow farn; + farn.pyrScale = pyrScale; + farn.polyN = polyN; + farn.polySigma = polySigma; + farn.flags = flags; + + cv::ocl::oclMat d_flowx, d_flowy; + cv::Mat flow, flowBuf, flowxBuf, flowyBuf; + + WARMUP_ON; + farn(d_frame0, d_frame1, d_flowx, d_flowy); + + if (useInitFlow) + { + cv::Mat flowxy[] = {cv::Mat(d_flowx), cv::Mat(d_flowy)}; + cv::merge(flowxy, 2, flow); + flow.copyTo(flowBuf); + flowxy[0].copyTo(flowxBuf); + flowxy[1].copyTo(flowyBuf); + + farn.flags |= cv::OPTFLOW_USE_INITIAL_FLOW; + farn(d_frame0, d_frame1, d_flowx, d_flowy); + } + WARMUP_OFF; + + cv::calcOpticalFlowFarneback( + frame0, frame1, flow, farn.pyrScale, farn.numLevels, farn.winSize, + farn.numIters, farn.polyN, farn.polySigma, farn.flags); + + std::vector flowxy; + cv::split(flow, flowxy); + + double diff0 = 0.0; + TestSystem::instance().setAccurate(ExceptedMatSimilar(flowxy[0], cv::Mat(d_flowx), 0.1, diff0)); + TestSystem::instance().setDiff(diff0); + double diff1 = 0.0; + TestSystem::instance().setAccurate(ExceptedMatSimilar(flowxy[1], cv::Mat(d_flowy), 0.1, diff1)); + TestSystem::instance().setDiff(diff1); + + if (useInitFlow) + { + cv::Mat flowx, flowy; + farn.flags = (flags | cv::OPTFLOW_USE_INITIAL_FLOW); + + CPU_ON; + cv::calcOpticalFlowFarneback( + frame0, frame1, flowBuf, farn.pyrScale, farn.numLevels, farn.winSize, + farn.numIters, farn.polyN, farn.polySigma, farn.flags); + CPU_OFF; + + GPU_ON; + farn(d_frame0, d_frame1, d_flowx, d_flowy); + GPU_OFF; + + GPU_FULL_ON; + d_frame0.upload(frame0); + d_frame1.upload(frame1); + d_flowx.upload(flowxBuf); + d_flowy.upload(flowyBuf); + farn(d_frame0, d_frame1, d_flowx, d_flowy); + d_flowx.download(flowx); + d_flowy.download(flowy); + GPU_FULL_OFF; + } + else + { + cv::Mat flow, flowx, flowy; + cv::ocl::oclMat d_flowx, d_flowy; + + farn.flags = flags; + + CPU_ON; + cv::calcOpticalFlowFarneback( + frame0, frame1, flow, farn.pyrScale, farn.numLevels, farn.winSize, + farn.numIters, farn.polyN, farn.polySigma, farn.flags); + CPU_OFF; + + GPU_ON; + farn(d_frame0, d_frame1, d_flowx, d_flowy); + GPU_OFF; + + GPU_FULL_ON; + d_frame0.upload(frame0); + d_frame1.upload(frame1); + farn(d_frame0, d_frame1, d_flowx, d_flowy); + d_flowx.download(flowx); + d_flowy.download(flowy); + GPU_FULL_OFF; + } + } + } + } +} + diff --git a/modules/ocl/src/opencl/optical_flow_farneback.cl b/modules/ocl/src/opencl/optical_flow_farneback.cl new file mode 100644 index 0000000000..2e5c6d9606 --- /dev/null +++ b/modules/ocl/src/opencl/optical_flow_farneback.cl @@ -0,0 +1,446 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2010-2012, Multicoreware, Inc., all rights reserved. +// Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved. +// Third party copyrights are property of their respective owners. +// +// @Authors +// Sen Liu, swjtuls1987@126.com +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistribution's of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistribution's in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other oclMaterials provided with the distribution. +// +// * The name of the copyright holders may not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// This software is provided by the copyright holders and contributors as is and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall the Intel Corporation or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + + +#define tx get_local_id(0) +#define ty get_local_id(1) +#define bx get_group_id(0) +#define bdx get_local_size(0) + +#define BORDER_SIZE 5 +#define MAX_KSIZE_HALF 100 + +#ifndef polyN +#define polyN 5 +#endif + +__kernel void polynomialExpansion(__global float * dst, + __global __const float * src, + __global __const float * c_g, + __global __const float * c_xg, + __global __const float * c_xxg, + __local float * smem, + const float4 ig, + const int height, const int width, + int dstStep, int srcStep) +{ + const int y = get_global_id(1); + const int x = bx * (bdx - 2*polyN) + tx - polyN; + + dstStep /= sizeof(*dst); + srcStep /= sizeof(*src); + + int xWarped; + __local float *row = smem + tx; + + if (y < height && y >= 0) + { + xWarped = min(max(x, 0), width - 1); + + row[0] = src[mad24(y, srcStep, xWarped)] * c_g[0]; + row[bdx] = 0.f; + row[2*bdx] = 0.f; + +#pragma unroll + for (int k = 1; k <= polyN; ++k) + { + float t0 = src[mad24(max(y - k, 0), srcStep, xWarped)]; + float t1 = src[mad24(min(y + k, height - 1), srcStep, xWarped)]; + + row[0] += c_g[k] * (t0 + t1); + row[bdx] += c_xg[k] * (t1 - t0); + row[2*bdx] += c_xxg[k] * (t0 + t1); + } + } + + barrier(CLK_LOCAL_MEM_FENCE); + + if (y < height && y >= 0 && tx >= polyN && tx + polyN < bdx && x < width) + { + float b1 = c_g[0] * row[0]; + float b3 = c_g[0] * row[bdx]; + float b5 = c_g[0] * row[2*bdx]; + float b2 = 0, b4 = 0, b6 = 0; + +#pragma unroll + for (int k = 1; k <= polyN; ++k) + { + b1 += (row[k] + row[-k]) * c_g[k]; + b4 += (row[k] + row[-k]) * c_xxg[k]; + b2 += (row[k] - row[-k]) * c_xg[k]; + b3 += (row[k + bdx] + row[-k + bdx]) * c_g[k]; + b6 += (row[k + bdx] - row[-k + bdx]) * c_xg[k]; + b5 += (row[k + 2*bdx] + row[-k + 2*bdx]) * c_g[k]; + } + + dst[mad24(y, dstStep, xWarped)] = b3*ig.s0; + dst[mad24(height + y, dstStep, xWarped)] = b2*ig.s0; + dst[mad24(2*height + y, dstStep, xWarped)] = b1*ig.s1 + b5*ig.s2; + dst[mad24(3*height + y, dstStep, xWarped)] = b1*ig.s1 + b4*ig.s2; + dst[mad24(4*height + y, dstStep, xWarped)] = b6*ig.s3; + } +} + +inline int idx_row_low(const int y, const int last_row) +{ + return abs(y) % (last_row + 1); +} + +inline int idx_row_high(const int y, const int last_row) +{ + return abs(last_row - abs(last_row - y)) % (last_row + 1); +} + +inline int idx_row(const int y, const int last_row) +{ + return idx_row_low(idx_row_high(y, last_row), last_row); +} + +inline int idx_col_low(const int x, const int last_col) +{ + return abs(x) % (last_col + 1); +} + +inline int idx_col_high(const int x, const int last_col) +{ + return abs(last_col - abs(last_col - x)) % (last_col + 1); +} + +inline int idx_col(const int x, const int last_col) +{ + return idx_col_low(idx_col_high(x, last_col), last_col); +} + +__kernel void gaussianBlur(__global float * dst, + __global const float * src, + __global const float * c_gKer, + __local float * smem, + const int height, const int width, + int dstStep, int srcStep, + const int ksizeHalf) +{ + const int y = get_global_id(1); + const int x = get_global_id(0); + + dstStep /= sizeof(*dst); + srcStep /= sizeof(*src); + + __local float *row = smem + ty * (bdx + 2*ksizeHalf); + + if (y < height) + { + // Vertical pass + for (int i = tx; i < bdx + 2*ksizeHalf; i += bdx) + { + int xExt = (int)(bx * bdx) + i - ksizeHalf; + xExt = idx_col(xExt, width - 1); + row[i] = src[mad24(y, srcStep, xExt)] * c_gKer[0]; + for (int j = 1; j <= ksizeHalf; ++j) + row[i] += (src[mad24(idx_row_low(y - j, height - 1), srcStep, xExt)] + + src[mad24(idx_row_high(y + j, height - 1), srcStep, xExt)]) * c_gKer[j]; + } + } + + barrier(CLK_LOCAL_MEM_FENCE); + + if (y < height && y >= 0 && x < width && x >= 0) + { + // Horizontal pass + row += tx + ksizeHalf; + float res = row[0] * c_gKer[0]; + for (int i = 1; i <= ksizeHalf; ++i) + res += (row[-i] + row[i]) * c_gKer[i]; + + dst[mad24(y, dstStep, x)] = res; + } +} + +__constant float c_border[BORDER_SIZE + 1] = { 0.14f, 0.14f, 0.4472f, 0.4472f, 0.4472f, 1.f }; + +__kernel void updateMatrices(__global float * M, + __global const float * flowx, __global const float * flowy, + __global const float * R0, __global const float * R1, + const int height, const int width, + int mStep, int xStep, int yStep, int R0Step, int R1Step) +{ + const int y = get_global_id(1); + const int x = get_global_id(0); + + mStep /= sizeof(*M); + xStep /= sizeof(*flowx); + yStep /= sizeof(*flowy); + R0Step /= sizeof(*R0); + R1Step /= sizeof(*R1); + + if (y < height && y >= 0 && x < width && x >= 0) + { + float dx = flowx[mad24(y, xStep, x)]; + float dy = flowy[mad24(y, yStep, x)]; + float fx = x + dx; + float fy = y + dy; + + int x1 = convert_int(floor(fx)); + int y1 = convert_int(floor(fy)); + fx -= x1; fy -= y1; + + float r2, r3, r4, r5, r6; + + if (x1 >= 0 && y1 >= 0 && x1 < width - 1 && y1 < height - 1) + { + float a00 = (1.f - fx) * (1.f - fy); + float a01 = fx * (1.f - fy); + float a10 = (1.f - fx) * fy; + float a11 = fx * fy; + + r2 = a00 * R1[mad24(y1, R1Step, x1)] + + a01 * R1[mad24(y1, R1Step, x1 + 1)] + + a10 * R1[mad24(y1 + 1, R1Step, x1)] + + a11 * R1[mad24(y1 + 1, R1Step, x1 + 1)]; + + r3 = a00 * R1[mad24(height + y1, R1Step, x1)] + + a01 * R1[mad24(height + y1, R1Step, x1 + 1)] + + a10 * R1[mad24(height + y1 + 1, R1Step, x1)] + + a11 * R1[mad24(height + y1 + 1, R1Step, x1 + 1)]; + + r4 = a00 * R1[mad24(2*height + y1, R1Step, x1)] + + a01 * R1[mad24(2*height + y1, R1Step, x1 + 1)] + + a10 * R1[mad24(2*height + y1 + 1, R1Step, x1)] + + a11 * R1[mad24(2*height + y1 + 1, R1Step, x1 + 1)]; + + r5 = a00 * R1[mad24(3*height + y1, R1Step, x1)] + + a01 * R1[mad24(3*height + y1, R1Step, x1 + 1)] + + a10 * R1[mad24(3*height + y1 + 1, R1Step, x1)] + + a11 * R1[mad24(3*height + y1 + 1, R1Step, x1 + 1)]; + + r6 = a00 * R1[mad24(4*height + y1, R1Step, x1)] + + a01 * R1[mad24(4*height + y1, R1Step, x1 + 1)] + + a10 * R1[mad24(4*height + y1 + 1, R1Step, x1)] + + a11 * R1[mad24(4*height + y1 + 1, R1Step, x1 + 1)]; + + r4 = (R0[mad24(2*height + y, R0Step, x)] + r4) * 0.5f; + r5 = (R0[mad24(3*height + y, R0Step, x)] + r5) * 0.5f; + r6 = (R0[mad24(4*height + y, R0Step, x)] + r6) * 0.25f; + } + else + { + r2 = r3 = 0.f; + r4 = R0[mad24(2*height + y, R0Step, x)]; + r5 = R0[mad24(3*height + y, R0Step, x)]; + r6 = R0[mad24(4*height + y, R0Step, x)] * 0.5f; + } + + r2 = (R0[mad24(y, R0Step, x)] - r2) * 0.5f; + r3 = (R0[mad24(height + y, R0Step, x)] - r3) * 0.5f; + + r2 += r4*dy + r6*dx; + r3 += r6*dy + r5*dx; + + float scale = + c_border[min(x, BORDER_SIZE)] * + c_border[min(y, BORDER_SIZE)] * + c_border[min(width - x - 1, BORDER_SIZE)] * + c_border[min(height - y - 1, BORDER_SIZE)]; + + r2 *= scale; r3 *= scale; r4 *= scale; + r5 *= scale; r6 *= scale; + + M[mad24(y, mStep, x)] = r4*r4 + r6*r6; + M[mad24(height + y, mStep, x)] = (r4 + r5)*r6; + M[mad24(2*height + y, mStep, x)] = r5*r5 + r6*r6; + M[mad24(3*height + y, mStep, x)] = r4*r2 + r6*r3; + M[mad24(4*height + y, mStep, x)] = r6*r2 + r5*r3; + } +} + +__kernel void boxFilter5(__global float * dst, + __global const float * src, + __local float * smem, + const int height, const int width, + int dstStep, int srcStep, + const int ksizeHalf) +{ + const int y = get_global_id(1); + const int x = get_global_id(0); + + const float boxAreaInv = 1.f / ((1 + 2*ksizeHalf) * (1 + 2*ksizeHalf)); + const int smw = bdx + 2*ksizeHalf; // shared memory "width" + __local float *row = smem + 5 * ty * smw; + + dstStep /= sizeof(*dst); + srcStep /= sizeof(*src); + + if (y < height) + { + // Vertical pass + for (int i = tx; i < bdx + 2*ksizeHalf; i += bdx) + { + int xExt = (int)(bx * bdx) + i - ksizeHalf; + xExt = min(max(xExt, 0), width - 1); + + #pragma unroll + for (int k = 0; k < 5; ++k) + row[k*smw + i] = src[mad24(k*height + y, srcStep, xExt)]; + + for (int j = 1; j <= ksizeHalf; ++j) + #pragma unroll + for (int k = 0; k < 5; ++k) + row[k*smw + i] += + src[mad24(k*height + max(y - j, 0), srcStep, xExt)] + + src[mad24(k*height + min(y + j, height - 1), srcStep, xExt)]; + } + } + + barrier(CLK_LOCAL_MEM_FENCE); + + if (y < height && y >= 0 && x < width && x >= 0) + { + // Horizontal pass + + row += tx + ksizeHalf; + float res[5]; + + #pragma unroll + for (int k = 0; k < 5; ++k) + res[k] = row[k*smw]; + + for (int i = 1; i <= ksizeHalf; ++i) + #pragma unroll + for (int k = 0; k < 5; ++k) + res[k] += row[k*smw - i] + row[k*smw + i]; + + #pragma unroll + for (int k = 0; k < 5; ++k) + dst[mad24(k*height + y, dstStep, x)] = res[k] * boxAreaInv; + } +} + +__kernel void updateFlow(__global float4 * flowx, __global float4 * flowy, + __global const float4 * M, + const int height, const int width, + int xStep, int yStep, int mStep) +{ + const int y = get_global_id(1); + const int x = get_global_id(0); + + xStep /= sizeof(*flowx); + yStep /= sizeof(*flowy); + mStep /= sizeof(*M); + + if (y < height && y >= 0 && x < width && x >= 0) + { + float4 g11 = M[mad24(y, mStep, x)]; + float4 g12 = M[mad24(height + y, mStep, x)]; + float4 g22 = M[mad24(2*height + y, mStep, x)]; + float4 h1 = M[mad24(3*height + y, mStep, x)]; + float4 h2 = M[mad24(4*height + y, mStep, x)]; + + float4 detInv = (float4)(1.f) / (g11*g22 - g12*g12 + (float4)(1e-3f)); + + flowx[mad24(y, xStep, x)] = (g11*h2 - g12*h1) * detInv; + flowy[mad24(y, yStep, x)] = (g22*h1 - g12*h2) * detInv; + } +} + +__kernel void gaussianBlur5(__global float * dst, + __global const float * src, + __global const float * c_gKer, + __local float * smem, + const int height, const int width, + int dstStep, int srcStep, + const int ksizeHalf) +{ + const int y = get_global_id(1); + const int x = get_global_id(0); + + const int smw = bdx + 2*ksizeHalf; // shared memory "width" + __local volatile float *row = smem + 5 * ty * smw; + + dstStep /= sizeof(*dst); + srcStep /= sizeof(*src); + + if (y < height) + { + // Vertical pass + for (int i = tx; i < bdx + 2*ksizeHalf; i += bdx) + { + int xExt = (int)(bx * bdx) + i - ksizeHalf; + xExt = idx_col(xExt, width - 1); + + #pragma unroll + for (int k = 0; k < 5; ++k) + row[k*smw + i] = src[mad24(k*height + y, srcStep, xExt)] * c_gKer[0]; + + for (int j = 1; j <= ksizeHalf; ++j) + #pragma unroll + for (int k = 0; k < 5; ++k) + row[k*smw + i] += + (src[mad24(k*height + idx_row_low(y - j, height - 1), srcStep, xExt)] + + src[mad24(k*height + idx_row_high(y + j, height - 1), srcStep, xExt)]) * c_gKer[j]; + } + } + + barrier(CLK_LOCAL_MEM_FENCE); + + if (y < height && y >= 0 && x < width && x >= 0) + { + // Horizontal pass + + row += tx + ksizeHalf; + float res[5]; + + #pragma unroll + for (int k = 0; k < 5; ++k) + res[k] = row[k*smw] * c_gKer[0]; + + for (int i = 1; i <= ksizeHalf; ++i) + #pragma unroll + for (int k = 0; k < 5; ++k) + res[k] += (row[k*smw - i] + row[k*smw + i]) * c_gKer[i]; + + #pragma unroll + for (int k = 0; k < 5; ++k) + dst[mad24(k*height + y, dstStep, x)] = res[k]; + } +} diff --git a/modules/ocl/src/optical_flow_farneback.cpp b/modules/ocl/src/optical_flow_farneback.cpp new file mode 100644 index 0000000000..6667eb7253 --- /dev/null +++ b/modules/ocl/src/optical_flow_farneback.cpp @@ -0,0 +1,507 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2010-2012, Multicoreware, Inc., all rights reserved. +// Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved. +// Third party copyrights are property of their respective owners. +// +// @Authors +// Sen Liu, swjtuls1987@126.com +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistribution's of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistribution's in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other oclMaterials provided with the distribution. +// +// * The name of the copyright holders may not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall the Intel Corporation or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + + +#include "precomp.hpp" +#include "opencv2/video/tracking.hpp" + +using namespace std; +using namespace cv; +using namespace cv::ocl; + +#define MIN_SIZE 32 + +namespace cv +{ + namespace ocl + { + ///////////////////////////OpenCL kernel strings/////////////////////////// + extern const char *optical_flow_farneback; + } +} + +namespace cv { namespace ocl { namespace optflow_farneback +{ + oclMat g; + oclMat xg; + oclMat xxg; + oclMat gKer; + + float ig[4]; + + inline int divUp(int total, int grain) + { + return (total + grain - 1) / grain; + } + + inline void setGaussianBlurKernel(const float *c_gKer, int ksizeHalf) + { + cv::Mat t_gKer(1, ksizeHalf + 1, CV_32FC1, const_cast(c_gKer)); + gKer.upload(t_gKer); + } + + void gaussianBlurOcl(const oclMat &src, int ksizeHalf, oclMat &dst) + { + string kernelName("gaussianBlur"); + size_t localThreads[3] = { 256, 1, 1 }; + size_t globalThreads[3] = { divUp(src.cols, localThreads[0]) * localThreads[0], src.rows, 1 }; + int smem_size = (localThreads[0] + 2*ksizeHalf) * sizeof(float); + + CV_Assert(dst.size() == src.size()); + std::vector< std::pair > args; + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&dst.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&src.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&gKer.data)); + args.push_back(std::make_pair(smem_size, (void *)NULL)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&dst.rows)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&dst.cols)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&dst.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&src.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&ksizeHalf)); + + openCLExecuteKernel(Context::getContext(), &optical_flow_farneback, kernelName, + globalThreads, localThreads, args, -1, -1); + } + + void polynomialExpansionOcl(const oclMat &src, int polyN, oclMat &dst) + { + string kernelName("polynomialExpansion"); + size_t localThreads[3] = { 256, 1, 1 }; + size_t globalThreads[3] = { divUp(src.cols, localThreads[0] - 2*polyN) * localThreads[0], src.rows, 1 }; + int smem_size = 3 * localThreads[0] * sizeof(float); + + std::vector< std::pair > args; + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&dst.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&src.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&g.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&xg.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&xxg.data)); + args.push_back(std::make_pair(smem_size, (void *)NULL)); + args.push_back(std::make_pair(sizeof(cl_float4), (void *)&ig)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&src.rows)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&src.cols)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&dst.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&src.step)); + + char opt [128]; + sprintf(opt, "-D polyN=%d", polyN); + + openCLExecuteKernel(Context::getContext(), &optical_flow_farneback, kernelName, + globalThreads, localThreads, args, -1, -1, opt); + } + + void updateMatricesOcl(const oclMat &flowx, const oclMat &flowy, const oclMat &R0, const oclMat &R1, oclMat &M) + { + string kernelName("updateMatrices"); + size_t localThreads[3] = { 32, 8, 1 }; + size_t globalThreads[3] = { divUp(flowx.cols, localThreads[0]) * localThreads[0], + divUp(flowx.rows, localThreads[1]) * localThreads[1], + 1 }; + + std::vector< std::pair > args; + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&M.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&flowx.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&flowy.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&R0.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&R1.data)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&flowx.rows)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&flowx.cols)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&M.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&flowx.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&flowy.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&R0.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&R1.step)); + + openCLExecuteKernel(Context::getContext(), &optical_flow_farneback, kernelName, + globalThreads, localThreads, args, -1, -1); + } + + void boxFilter5Ocl(const oclMat &src, int ksizeHalf, oclMat &dst) + { + string kernelName("boxFilter5"); + int height = src.rows / 5; + size_t localThreads[3] = { 256, 1, 1 }; + size_t globalThreads[3] = { divUp(src.cols, localThreads[0]) * localThreads[0], height, 1 }; + int smem_size = (localThreads[0] + 2*ksizeHalf) * 5 * sizeof(float); + + std::vector< std::pair > args; + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&dst.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&src.data)); + args.push_back(std::make_pair(smem_size, (void *)NULL)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&height)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&src.cols)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&dst.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&src.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&ksizeHalf)); + + openCLExecuteKernel(Context::getContext(), &optical_flow_farneback, kernelName, + globalThreads, localThreads, args, -1, -1); + } + + void updateFlowOcl(const oclMat &M, oclMat &flowx, oclMat &flowy) + { + string kernelName("updateFlow"); + int cols = divUp(flowx.cols, 4); + size_t localThreads[3] = { 32, 8, 1 }; + size_t globalThreads[3] = { divUp(cols, localThreads[0]) * localThreads[0], + divUp(flowx.rows, localThreads[1]) * localThreads[0], + 1 }; + + std::vector< std::pair > args; + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&flowx.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&flowy.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&M.data)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&flowx.rows)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&cols)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&flowx.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&flowy.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&M.step)); + + openCLExecuteKernel(Context::getContext(), &optical_flow_farneback, kernelName, + globalThreads, localThreads, args, -1, -1); + } + + void gaussianBlur5Ocl(const oclMat &src, int ksizeHalf, oclMat &dst) + { + string kernelName("gaussianBlur5"); + int height = src.rows / 5; + int width = src.cols; + size_t localThreads[3] = { 256, 1, 1 }; + size_t globalThreads[3] = { divUp(width, localThreads[0]) * localThreads[0], height, 1 }; + int smem_size = (localThreads[0] + 2*ksizeHalf) * 5 * sizeof(float); + + std::vector< std::pair > args; + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&dst.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&src.data)); + args.push_back(std::make_pair(sizeof(cl_mem), (void *)&gKer.data)); + args.push_back(std::make_pair(smem_size, (void *)NULL)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&height)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&width)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&dst.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&src.step)); + args.push_back(std::make_pair(sizeof(cl_int), (void *)&ksizeHalf)); + + openCLExecuteKernel(Context::getContext(), &optical_flow_farneback, kernelName, + globalThreads, localThreads, args, -1, -1); + } +}}} // namespace cv { namespace ocl { namespace optflow_farneback + +static oclMat allocMatFromBuf(int rows, int cols, int type, oclMat &mat) +{ + if (!mat.empty() && mat.type() == type && mat.rows >= rows && mat.cols >= cols) + return mat(Rect(0, 0, cols, rows)); + return mat = oclMat(rows, cols, type); +} + +void cv::ocl::FarnebackOpticalFlow::prepareGaussian( + int n, double sigma, float *g, float *xg, float *xxg, + double &ig11, double &ig03, double &ig33, double &ig55) +{ + double s = 0.; + for (int x = -n; x <= n; x++) + { + g[x] = (float)std::exp(-x*x/(2*sigma*sigma)); + s += g[x]; + } + + s = 1./s; + for (int x = -n; x <= n; x++) + { + g[x] = (float)(g[x]*s); + xg[x] = (float)(x*g[x]); + xxg[x] = (float)(x*x*g[x]); + } + + Mat_ G(6, 6); + G.setTo(0); + + for (int y = -n; y <= n; y++) + { + for (int x = -n; x <= n; x++) + { + G(0,0) += g[y]*g[x]; + G(1,1) += g[y]*g[x]*x*x; + G(3,3) += g[y]*g[x]*x*x*x*x; + G(5,5) += g[y]*g[x]*x*x*y*y; + } + } + + //G[0][0] = 1.; + G(2,2) = G(0,3) = G(0,4) = G(3,0) = G(4,0) = G(1,1); + G(4,4) = G(3,3); + G(3,4) = G(4,3) = G(5,5); + + // invG: + // [ x e e ] + // [ y ] + // [ y ] + // [ e z ] + // [ e z ] + // [ u ] + Mat_ invG = G.inv(DECOMP_CHOLESKY); + + ig11 = invG(1,1); + ig03 = invG(0,3); + ig33 = invG(3,3); + ig55 = invG(5,5); +} + +void cv::ocl::FarnebackOpticalFlow::setPolynomialExpansionConsts(int n, double sigma) +{ + vector buf(n*6 + 3); + float* g = &buf[0] + n; + float* xg = g + n*2 + 1; + float* xxg = xg + n*2 + 1; + + if (sigma < FLT_EPSILON) + sigma = n*0.3; + + double ig11, ig03, ig33, ig55; + prepareGaussian(n, sigma, g, xg, xxg, ig11, ig03, ig33, ig55); + + cv::Mat t_g(1, n + 1, CV_32FC1, g); + cv::Mat t_xg(1, n + 1, CV_32FC1, xg); + cv::Mat t_xxg(1, n + 1, CV_32FC1, xxg); + + optflow_farneback::g.upload(t_g); + optflow_farneback::xg.upload(t_xg); + optflow_farneback::xxg.upload(t_xxg); + + optflow_farneback::ig[0] = static_cast(ig11); + optflow_farneback::ig[1] = static_cast(ig03); + optflow_farneback::ig[2] = static_cast(ig33); + optflow_farneback::ig[3] = static_cast(ig55); +} + +void cv::ocl::FarnebackOpticalFlow::updateFlow_boxFilter( + const oclMat& R0, const oclMat& R1, oclMat& flowx, oclMat &flowy, + oclMat& M, oclMat &bufM, int blockSize, bool updateMatrices) +{ + optflow_farneback::boxFilter5Ocl(M, blockSize/2, bufM); + + swap(M, bufM); + + finish(); + + optflow_farneback::updateFlowOcl(M, flowx, flowy); + + if (updateMatrices) + optflow_farneback::updateMatricesOcl(flowx, flowy, R0, R1, M); +} + + +void cv::ocl::FarnebackOpticalFlow::updateFlow_gaussianBlur( + const oclMat& R0, const oclMat& R1, oclMat& flowx, oclMat& flowy, + oclMat& M, oclMat &bufM, int blockSize, bool updateMatrices) +{ + optflow_farneback::gaussianBlur5Ocl(M, blockSize/2, bufM); + + swap(M, bufM); + + optflow_farneback::updateFlowOcl(M, flowx, flowy); + + if (updateMatrices) + optflow_farneback::updateMatricesOcl(flowx, flowy, R0, R1, M); +} + + +void cv::ocl::FarnebackOpticalFlow::operator ()( + const oclMat &frame0, const oclMat &frame1, oclMat &flowx, oclMat &flowy) +{ + CV_Assert(frame0.channels() == 1 && frame1.channels() == 1); + CV_Assert(frame0.size() == frame1.size()); + CV_Assert(polyN == 5 || polyN == 7); + CV_Assert(!fastPyramids || std::abs(pyrScale - 0.5) < 1e-6); + + Size size = frame0.size(); + oclMat prevFlowX, prevFlowY, curFlowX, curFlowY; + + flowx.create(size, CV_32F); + flowy.create(size, CV_32F); + oclMat flowx0 = flowx; + oclMat flowy0 = flowy; + + // Crop unnecessary levels + double scale = 1; + int numLevelsCropped = 0; + for (; numLevelsCropped < numLevels; numLevelsCropped++) + { + scale *= pyrScale; + if (size.width*scale < MIN_SIZE || size.height*scale < MIN_SIZE) + break; + } + + frame0.convertTo(frames_[0], CV_32F); + frame1.convertTo(frames_[1], CV_32F); + + if (fastPyramids) + { + // Build Gaussian pyramids using pyrDown() + pyramid0_.resize(numLevelsCropped + 1); + pyramid1_.resize(numLevelsCropped + 1); + pyramid0_[0] = frames_[0]; + pyramid1_[0] = frames_[1]; + for (int i = 1; i <= numLevelsCropped; ++i) + { + pyrDown(pyramid0_[i - 1], pyramid0_[i]); + pyrDown(pyramid1_[i - 1], pyramid1_[i]); + } + } + + setPolynomialExpansionConsts(polyN, polySigma); + + for (int k = numLevelsCropped; k >= 0; k--) + { + scale = 1; + for (int i = 0; i < k; i++) + scale *= pyrScale; + + double sigma = (1./scale - 1) * 0.5; + int smoothSize = cvRound(sigma*5) | 1; + smoothSize = std::max(smoothSize, 3); + + int width = cvRound(size.width*scale); + int height = cvRound(size.height*scale); + + if (fastPyramids) + { + width = pyramid0_[k].cols; + height = pyramid0_[k].rows; + } + + if (k > 0) + { + curFlowX.create(height, width, CV_32F); + curFlowY.create(height, width, CV_32F); + } + else + { + curFlowX = flowx0; + curFlowY = flowy0; + } + + if (!prevFlowX.data) + { + if (flags & cv::OPTFLOW_USE_INITIAL_FLOW) + { + resize(flowx0, curFlowX, Size(width, height), 0, 0, INTER_LINEAR); + resize(flowy0, curFlowY, Size(width, height), 0, 0, INTER_LINEAR); + multiply(scale, curFlowX, curFlowX); + multiply(scale, curFlowY, curFlowY); + } + else + { + curFlowX.setTo(0); + curFlowY.setTo(0); + } + } + else + { + resize(prevFlowX, curFlowX, Size(width, height), 0, 0, INTER_LINEAR); + resize(prevFlowY, curFlowY, Size(width, height), 0, 0, INTER_LINEAR); + multiply(1./pyrScale, curFlowX, curFlowX); + multiply(1./pyrScale, curFlowY, curFlowY); + } + + oclMat M = allocMatFromBuf(5*height, width, CV_32F, M_); + oclMat bufM = allocMatFromBuf(5*height, width, CV_32F, bufM_); + oclMat R[2] = + { + allocMatFromBuf(5*height, width, CV_32F, R_[0]), + allocMatFromBuf(5*height, width, CV_32F, R_[1]) + }; + + if (fastPyramids) + { + optflow_farneback::polynomialExpansionOcl(pyramid0_[k], polyN, R[0]); + optflow_farneback::polynomialExpansionOcl(pyramid1_[k], polyN, R[1]); + } + else + { + oclMat blurredFrame[2] = + { + allocMatFromBuf(size.height, size.width, CV_32F, blurredFrame_[0]), + allocMatFromBuf(size.height, size.width, CV_32F, blurredFrame_[1]) + }; + oclMat pyrLevel[2] = + { + allocMatFromBuf(height, width, CV_32F, pyrLevel_[0]), + allocMatFromBuf(height, width, CV_32F, pyrLevel_[1]) + }; + + Mat g = getGaussianKernel(smoothSize, sigma, CV_32F); + optflow_farneback::setGaussianBlurKernel(g.ptr(smoothSize/2), smoothSize/2); + + for (int i = 0; i < 2; i++) + { + optflow_farneback::gaussianBlurOcl(frames_[i], smoothSize/2, blurredFrame[i]); + resize(blurredFrame[i], pyrLevel[i], Size(width, height), INTER_LINEAR); + optflow_farneback::polynomialExpansionOcl(pyrLevel[i], polyN, R[i]); + } + } + + optflow_farneback::updateMatricesOcl(curFlowX, curFlowY, R[0], R[1], M); + + if (flags & OPTFLOW_FARNEBACK_GAUSSIAN) + { + Mat g = getGaussianKernel(winSize, winSize/2*0.3f, CV_32F); + optflow_farneback::setGaussianBlurKernel(g.ptr(winSize/2), winSize/2); + } + for (int i = 0; i < numIters; i++) + { + if (flags & OPTFLOW_FARNEBACK_GAUSSIAN) + updateFlow_gaussianBlur(R[0], R[1], curFlowX, curFlowY, M, bufM, winSize, i < numIters-1); + else + updateFlow_boxFilter(R[0], R[1], curFlowX, curFlowY, M, bufM, winSize, i < numIters-1); + } + + prevFlowX = curFlowX; + prevFlowY = curFlowY; + } + + flowx = curFlowX; + flowy = curFlowY; +} + diff --git a/modules/ocl/test/test_optflow.cpp b/modules/ocl/test/test_optflow.cpp index 0121be8f9e..34adb352c2 100644 --- a/modules/ocl/test/test_optflow.cpp +++ b/modules/ocl/test/test_optflow.cpp @@ -272,6 +272,78 @@ TEST_P(Sparse, Mat) INSTANTIATE_TEST_CASE_P(OCL_Video, Sparse, Combine( Values(false, true), Values(false, true))); +////////////////////////////////////////////////////// +// FarnebackOpticalFlow + +namespace +{ + IMPLEMENT_PARAM_CLASS(PyrScale, double) + IMPLEMENT_PARAM_CLASS(PolyN, int) + CV_FLAGS(FarnebackOptFlowFlags, 0, OPTFLOW_FARNEBACK_GAUSSIAN) + IMPLEMENT_PARAM_CLASS(UseInitFlow, bool) +} + +PARAM_TEST_CASE(Farneback, PyrScale, PolyN, FarnebackOptFlowFlags, UseInitFlow) +{ + double pyrScale; + int polyN; + int flags; + bool useInitFlow; + + virtual void SetUp() + { + pyrScale = GET_PARAM(0); + polyN = GET_PARAM(1); + flags = GET_PARAM(2); + useInitFlow = GET_PARAM(3); + } +}; + +TEST_P(Farneback, Accuracy) +{ + cv::Mat frame0 = imread(workdir + "/rubberwhale1.png", cv::IMREAD_GRAYSCALE); + ASSERT_FALSE(frame0.empty()); + + cv::Mat frame1 = imread(workdir + "/rubberwhale2.png", cv::IMREAD_GRAYSCALE); + ASSERT_FALSE(frame1.empty()); + + double polySigma = polyN <= 5 ? 1.1 : 1.5; + + cv::ocl::FarnebackOpticalFlow farn; + farn.pyrScale = pyrScale; + farn.polyN = polyN; + farn.polySigma = polySigma; + farn.flags = flags; + + cv::ocl::oclMat d_flowx, d_flowy; + farn(oclMat(frame0), oclMat(frame1), d_flowx, d_flowy); + + cv::Mat flow; + if (useInitFlow) + { + cv::Mat flowxy[] = {cv::Mat(d_flowx), cv::Mat(d_flowy)}; + cv::merge(flowxy, 2, flow); + + farn.flags |= cv::OPTFLOW_USE_INITIAL_FLOW; + farn(oclMat(frame0), oclMat(frame1), d_flowx, d_flowy); + } + + cv::calcOpticalFlowFarneback( + frame0, frame1, flow, farn.pyrScale, farn.numLevels, farn.winSize, + farn.numIters, farn.polyN, farn.polySigma, farn.flags); + + std::vector flowxy; + cv::split(flow, flowxy); + + EXPECT_MAT_SIMILAR(flowxy[0], d_flowx, 0.1); + EXPECT_MAT_SIMILAR(flowxy[1], d_flowy, 0.1); +} + +INSTANTIATE_TEST_CASE_P(OCL_Video, Farneback, testing::Combine( + testing::Values(PyrScale(0.3), PyrScale(0.5), PyrScale(0.8)), + testing::Values(PolyN(5), PolyN(7)), + testing::Values(FarnebackOptFlowFlags(0), FarnebackOptFlowFlags(cv::OPTFLOW_FARNEBACK_GAUSSIAN)), + testing::Values(UseInitFlow(false), UseInitFlow(true)))); #endif // HAVE_OPENCL