diff --git a/modules/photo/src/denoising.cpp b/modules/photo/src/denoising.cpp index c9a711b325..a3aa0f41f9 100644 --- a/modules/photo/src/denoising.cpp +++ b/modules/photo/src/denoising.cpp @@ -40,14 +40,17 @@ //M*/ #include "precomp.hpp" -#include "opencv2/photo.hpp" -#include "opencv2/imgproc.hpp" + #include "fast_nlmeans_denoising_invoker.hpp" #include "fast_nlmeans_multi_denoising_invoker.hpp" +#include "fast_nlmeans_denoising_opencl.hpp" void cv::fastNlMeansDenoising( InputArray _src, OutputArray _dst, float h, int templateWindowSize, int searchWindowSize) { + CV_OCL_RUN(_src.dims() <= 2 && (_src.isUMat() || _dst.isUMat()), + ocl_fastNlMeansDenoising(_src, _dst, h, templateWindowSize, searchWindowSize)) + Mat src = _src.getMat(); _dst.create(src.size(), src.type()); Mat dst = _dst.getMat(); diff --git a/modules/photo/src/fast_nlmeans_denoising_opencl.hpp b/modules/photo/src/fast_nlmeans_denoising_opencl.hpp new file mode 100644 index 0000000000..404a59ee09 --- /dev/null +++ b/modules/photo/src/fast_nlmeans_denoising_opencl.hpp @@ -0,0 +1,128 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. + +// Copyright (C) 2014, Advanced Micro Devices, Inc., all rights reserved. +// Third party copyrights are property of their respective owners. + +#ifndef __OPENCV_FAST_NLMEANS_DENOISING_OPENCL_HPP__ +#define __OPENCV_FAST_NLMEANS_DENOISING_OPENCL_HPP__ + +#include "precomp.hpp" + +#define CV_OPENCL_RUN_ASSERT +#include "opencl_kernels.hpp" + +namespace cv { + +enum +{ + BLOCK_ROWS = 32, + BLOCK_COLS = 128, + CTA_SIZE = 128 +}; + +static inline int getNearestPowerOf2(int value) +{ + int p = 0; + while (1 << p < value) + ++p; + return p; +} + +static int divUp(int a, int b) +{ + return (a + b - 1) / b; +} + +static bool ocl_calcAlmostDist2Weight(UMat & almostDist2Weight, int searchWindowSize, int templateWindowSize, float h, int cn, + int & almostTemplateWindowSizeSqBinShift) +{ + const int maxEstimateSumValue = searchWindowSize * searchWindowSize * 255; + int fixedPointMult = std::numeric_limits::max() / maxEstimateSumValue; + + // precalc weight for every possible l2 dist between blocks + // additional optimization of precalced weights to replace division(averaging) by binary shift + CV_Assert(templateWindowSize <= 46340); // sqrt(INT_MAX) + int templateWindowSizeSq = templateWindowSize * templateWindowSize; + almostTemplateWindowSizeSqBinShift = getNearestPowerOf2(templateWindowSizeSq); + float almostDist2ActualDistMultiplier = (float)(1 << almostTemplateWindowSizeSqBinShift) / templateWindowSizeSq; + + const float WEIGHT_THRESHOLD = 1e-3f; + int maxDist = 255 * 255 * cn; + int almostMaxDist = (int)(maxDist / almostDist2ActualDistMultiplier + 1); + float den = 1.0f / (h * h * cn); + + almostDist2Weight.create(1, almostMaxDist, CV_32SC1); + + ocl::Kernel k("calcAlmostDist2Weight", ocl::photo::nlmeans_oclsrc, + "-D OP_CALC_WEIGHTS"); + if (k.empty()) + return false; + + k.args(ocl::KernelArg::PtrWriteOnly(almostDist2Weight), almostMaxDist, + almostDist2ActualDistMultiplier, fixedPointMult, den, WEIGHT_THRESHOLD); + + size_t globalsize[1] = { almostMaxDist }; + return k.run(1, globalsize, NULL, false); +} + +static bool ocl_fastNlMeansDenoising(InputArray _src, OutputArray _dst, float h, + int templateWindowSize, int searchWindowSize) +{ + int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type); + Size size = _src.size(); + + if ( !(depth == CV_8U && cn <= 4 && cn != 3) ) + return false; + + int templateWindowHalfWize = templateWindowSize / 2; + int searchWindowHalfSize = searchWindowSize / 2; + templateWindowSize = templateWindowHalfWize * 2 + 1; + searchWindowSize = searchWindowHalfSize * 2 + 1; + int nblocksx = divUp(size.width, BLOCK_COLS), nblocksy = divUp(size.height, BLOCK_ROWS); + int almostTemplateWindowSizeSqBinShift = -1; + + char cvt[2][40]; + String opts = format("-D OP_CALC_FASTNLMEANS -D TEMPLATE_SIZE=%d -D SEARCH_SIZE=%d" + " -D uchar_t=%s -D int_t=%s -D BLOCK_COLS=%d -D BLOCK_ROWS=%d" + " -D CTA_SIZE=%d -D TEMPLATE_SIZE2=%d -D SEARCH_SIZE2=%d" + " -D convert_int_t=%s -D cn=%d -D CTA_SIZE2=%d -D convert_uchar_t=%s", + templateWindowSize, searchWindowSize, ocl::typeToStr(type), + ocl::typeToStr(CV_32SC(cn)), BLOCK_COLS, BLOCK_ROWS, CTA_SIZE, + templateWindowHalfWize, searchWindowHalfSize, + ocl::convertTypeStr(CV_8U, CV_32S, cn, cvt[0]), cn, + CTA_SIZE >> 1, ocl::convertTypeStr(CV_32S, CV_8U, cn, cvt[1])); + + ocl::Kernel k("fastNlMeansDenoising", ocl::photo::nlmeans_oclsrc, opts); + if (k.empty()) + return false; + + UMat almostDist2Weight; + if (!ocl_calcAlmostDist2Weight(almostDist2Weight, searchWindowSize, templateWindowSize, h, cn, + almostTemplateWindowSizeSqBinShift)) + return false; + CV_Assert(almostTemplateWindowSizeSqBinShift >= 0); + + UMat srcex; + int borderSize = searchWindowHalfSize + templateWindowHalfWize; + copyMakeBorder(_src, srcex, borderSize, borderSize, borderSize, borderSize, BORDER_DEFAULT); + + _dst.create(size, type); + UMat dst = _dst.getUMat(); + + Size upColSumSize(size.width, searchWindowSize * searchWindowSize * nblocksy); + Size colSumSize(nblocksx * templateWindowSize, searchWindowSize * searchWindowSize * nblocksy); + UMat buffer(upColSumSize + colSumSize, CV_32SC(cn)); + + k.args(ocl::KernelArg::ReadOnlyNoSize(srcex), ocl::KernelArg::WriteOnly(dst), + ocl::KernelArg::PtrReadOnly(almostDist2Weight), nblocksy, nblocksx, + ocl::KernelArg::PtrReadOnly(buffer), almostTemplateWindowSizeSqBinShift); + + size_t globalsize[2] = { nblocksx, nblocksy }, localsize[2] = { CTA_SIZE, 1 }; + return k.run(2, globalsize, localsize, false); +} + +} + +#endif diff --git a/modules/photo/src/opencl/nlmeans.cl b/modules/photo/src/opencl/nlmeans.cl new file mode 100644 index 0000000000..1b40517abe --- /dev/null +++ b/modules/photo/src/opencl/nlmeans.cl @@ -0,0 +1,249 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. + +// Copyright (C) 2014, Advanced Micro Devices, Inc., all rights reserved. +// Third party copyrights are property of their respective owners. + +#ifdef OP_CALC_WEIGHTS + +__kernel void calcAlmostDist2Weight(__global int * almostDist2Weight, int almostMaxDist, + float almostDist2ActualDistMultiplier, int fixedPointMult, + float den, float WEIGHT_THRESHOLD) +{ + int almostDist = get_global_id(0); + + if (almostDist < almostMaxDist) + { + float dist = almostDist * almostDist2ActualDistMultiplier; + int weight = convert_int_sat_rte(fixedPointMult * exp(-dist * den)); + + if (weight < WEIGHT_THRESHOLD * fixedPointMult) + weight = 0; + + almostDist2Weight[almostDist] = weight; + } +} + +#elif defined OP_CALC_FASTNLMEANS + +#define SEARCH_SIZE_SQ (SEARCH_SIZE * SEARCH_SIZE) + +inline int_t calcDist(uchar_t a, uchar_t b) +{ + int_t diff = convert_int_t(a) -convert_int_t(b); + return diff * diff; +} + +inline void calcFirstElementInRow(__global const uchar * src, int src_step, int src_offset, + __local int_t * dists, int y, int x, int id, + __global int_t * col_dists, __global int_t * up_col_dists) +{ + int sx = x - SEARCH_SIZE2, sy = y - SEARCH_SIZE2; + + for (int i = 0, size = SEARCH_SIZE_SQ; i < size; i += CTA_SIZE) + { + int_t dist = (int_t)(0), value; + + sx += i % SEARCH_SIZE; + sy += i / SEARCH_SIZE; + + __global const uchar_t * src_template = (__global const uchar_t *)(src + mad24(sy, src_step, mad24(cn, x, src_offset))); + __global const uchar_t * src_current = (__global const uchar_t *)(src + mad24(y, src_step, mad24(cn, x, src_offset))); + __global int_t * col_dists_current = col_dists + i * TEMPLATE_SIZE; + + #pragma unroll + for (int j = 0; j < TEMPLATE_SIZE; ++j) + col_dists_current[j] = (int_t)(0); + + #pragma unroll + for (int ty = -TEMPLATE_SIZE2; ty <= TEMPLATE_SIZE2; ++ty) + { + #pragma unroll + for (int tx = -TEMPLATE_SIZE2; tx <= TEMPLATE_SIZE2; ++tx) + { + value = calcDist(src_template[tx], src_current[tx]); + + col_dists_current[tx + TEMPLATE_SIZE2] += value; + dist += value; + } + + src_current += src_step; + src_template += src_step; + } + + dists[i] = dist; + up_col_dists[i] = col_dists[TEMPLATE_SIZE - 1]; + } +} + +inline void calcElementInFirstRow(__global const uchar * src, int src_step, int src_offset, + __local int_t * dists, int y, int x, int id, int first, + __global int_t * col_dists, __global int_t * up_col_dists) +{ + x += TEMPLATE_SIZE2; + int sx = x - SEARCH_SIZE2, sy = y - SEARCH_SIZE2; + + for (int i = 0, size = SEARCH_SIZE_SQ; i < size; i += CTA_SIZE) + { + sx += i % SEARCH_SIZE; + sy += i / SEARCH_SIZE; + + __global const uchar_t * src_current = (__global const uchar_t *)(src + mad24(y, src_step, mad24(cn, x, src_offset))); + __global const uchar_t * src_template = (__global const uchar_t *)(src + mad24(sy, src_step, mad24(cn, x, src_offset))); + __global int_t * col_dists_current = col_dists + TEMPLATE_SIZE * i; + + int_t value; + dists[id] -= col_dists_current[first]; + col_dists_current[first] = (int_t)(0); + + #pragma unroll + for (int ty = -TEMPLATE_SIZE2; ty <= TEMPLATE_SIZE2; ++ty) + { + value = calcDist(src_current[0], src_template[0]); + col_dists_current[first] += value; + + src_current += src_step; + src_template += src_step; + } + + dists[id] += col_dists_current[first]; + up_col_dists[id] = col_dists_current[first]; + } +} + +inline void calcElement(__global const uchar * src, int src_step, int src_offset, + __local int_t * dists, int y, int x, int id, int first, + __global int_t * col_dists, __global int_t * up_col_dists) +{ + int sx_up = x + TEMPLATE_SIZE2, sy_up = y - TEMPLATE_SIZE2 - 1; + int sx_down = x + TEMPLATE_SIZE2, sy_down = y + TEMPLATE_SIZE2; + + uchar_t up_value = *(__global const uchar_t *)(src + mad24(sy_up, src_step, mad24(cn, sx_up, src_offset))); + uchar_t down_value = *(__global const uchar_t *)(src + mad24(sy_down, src_step, mad24(cn, sx_down, src_offset))); + + for (int i = 0, size = SEARCH_SIZE_SQ; i < size; i += CTA_SIZE) + { + int wx = i % SEARCH_SIZE; + int wy = i / SEARCH_SIZE; + + sx_up += wx, sx_down += wx; + sy_up += wy, sy_down += wy; + + uchar_t up_value_t = *(__global const uchar_t *)(src + mad24(sy_up, src_step, mad24(cn, sx_up, src_offset))); + uchar_t down_value_t = *(__global const uchar_t *)(src + mad24(sy_down, src_step, mad24(cn, sx_down, src_offset))); + + __global int_t * col_dists_current = col_dists + i * TEMPLATE_SIZE; + __global int_t * up_col_dists_current = up_col_dists + i; + + dists[i] -= col_dists_current[first]; + col_dists_current[first] = up_col_dists_current[id] + calcDist(down_value, down_value_t) - calcDist(up_value, up_value_t); + dists[i] += col_dists_current[first]; + up_col_dists_current[id] = col_dists_current[first]; + } +} + +inline void convolveWindow(__global const uchar * src, int src_step, int src_offset, + __local int * dists, __global const int * almostDist2Weight, + __global uchar * dst, int dst_step, int dst_offset, + int y, int x, int id, __local int * weights_local, + __local int * weighted_sum_local, int almostTemplateWindowSizeSqBinShift) +{ + int sx = x - SEARCH_SIZE2, sy = y - SEARCH_SIZE2, weights = 0; + int_t weighted_sum = (int_t)(0); + + for (int i = 0, size = SEARCH_SIZE_SQ; i < size; i += id) + { + int src_index = mad24(sy + i / SEARCH_SIZE, src_step, (i % SEARCH_SIZE + sx) * cn + src_offset); + __global const uchar_t * src_search = (__global const uchar_t *)(src + src_index); + + int almostAvgDist = dists[i] >> almostTemplateWindowSizeSqBinShift; + int weight = almostDist2Weight[almostAvgDist]; + + weights += weight; + weighted_sum += (int_t)(weight) * convert_int_t(src_search[0]); + } + + if (id >= CTA_SIZE2) + { + weights_local[id - CTA_SIZE2] = weights; + weighted_sum_local[id - CTA_SIZE2] = weighted_sum; + } + barrier(CLK_LOCAL_MEM_FENCE); + + if (id < CTA_SIZE2) + { + weights_local[id] += weights; + weighted_sum_local[id] += weighted_sum; + } + barrier(CLK_LOCAL_MEM_FENCE); + + for (int lsize = CTA_SIZE2 >> 1; lsize >= 4; lsize >>= 1) + { + if (id < lsize) + { + int id2 = lsize + id; + weights_local[id] = weights + weights_local[id2]; + weighted_sum_local[id] = weighted_sum + weighted_sum_local[id2]; + } + barrier(CLK_LOCAL_MEM_FENCE); + } + + if (id == 0) + { + int dst_index = mad24(y, dst_step, dst_offset + x * cn); + + int_t weights_local_0 = (int_t)(weights_local[0] + weights_local[1] + weights_local[2] + weights_local[3]); + int_t weighted_sum_local_0 = weighted_sum_local[0] + weighted_sum_local[1] + weighted_sum_local[2] + weighted_sum_local[3]; + + *(__global uchar_t *)(dst + dst_index) = convert_uchar_t((weighted_sum_local_0 + weights_local_0 >> 1) / weights_local_0); + } +} + +__kernel void fastNlMeansDenoising(__global const uchar * src, int src_step, int src_offset, + __global uchar * dst, int dst_step, int dst_offset, int dst_rows, int dst_cols, + __global const int * almostDist2Weight, int nblocksy, int nblocksx, + __global uchar * buffer, int almostTemplateWindowSizeSqBinShift) +{ + int block_x = get_global_id(0); + int block_y = get_global_id(1); + int id = get_local_id(0), first; + + __local int_t dists[SEARCH_SIZE_SQ], weighted_sum[CTA_SIZE2]; + __local int weights[CTA_SIZE2]; + + int block_data_start = mad24(block_y, nblocksx, block_x) * SEARCH_SIZE_SQ * (TEMPLATE_SIZE + BLOCK_COLS); + __global int_t * col_dists = (__global int_t *)(buffer + block_data_start * sizeof(int_t)); + __global int_t * up_col_dists = (__global int_t *)(buffer + sizeof(int_t) * (block_data_start + SEARCH_SIZE_SQ * TEMPLATE_SIZE)); + + if (block_x < nblocksx && block_y < nblocksy) + { + int x0 = block_x * BLOCK_COLS, x1 = min(x0 + BLOCK_COLS, dst_cols); + int y0 = block_y * BLOCK_ROWS, y1 = min(y0 + BLOCK_ROWS, dst_rows); + + for (int y = y0; y < y1; ++y) + for (int x = x0; x < x1; ++x) + { + if (x == x0) + { + calcFirstElementInRow(src, src_step, src_offset, dists, y, x, id, col_dists, up_col_dists); + first = 0; + } + else + { + if (y == y0) + calcElementInFirstRow(src, src_step, src_offset, dists, y, x, id, first, col_dists, up_col_dists); + else + { + calcElement(src, src_step, src_offset, dists, y, x, id, first, col_dists, up_col_dists); + first = (first + 1) % TEMPLATE_SIZE; + } + + convolveWindow(src, src_step, src_offset, dists, almostDist2Weight, dst, dst_step, dst_offset, + y, x, id, weights, weighted_sum, almostTemplateWindowSizeSqBinShift); + } + } + } +} + +#endif diff --git a/modules/photo/src/precomp.hpp b/modules/photo/src/precomp.hpp index 38ac3ffcda..4355f13c00 100644 --- a/modules/photo/src/precomp.hpp +++ b/modules/photo/src/precomp.hpp @@ -46,6 +46,8 @@ #include "opencv2/core/private.hpp" #include "opencv2/core/utility.hpp" #include "opencv2/photo.hpp" +#include "opencv2/core/ocl.hpp" +#include "opencv2/imgproc.hpp" #ifdef HAVE_TEGRA_OPTIMIZATION #include "opencv2/photo/photo_tegra.hpp" diff --git a/modules/photo/test/ocl/test_denoising.cpp b/modules/photo/test/ocl/test_denoising.cpp new file mode 100644 index 0000000000..0554b30670 --- /dev/null +++ b/modules/photo/test/ocl/test_denoising.cpp @@ -0,0 +1,105 @@ +/*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) 2000-2008, Intel Corporation, all rights reserved. +// Copyright (C) 2009, Willow Garage Inc., all rights reserved. +// Third party copyrights are property of their respective owners. +// +// 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 materials 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 "test_precomp.hpp" +#include "opencv2/ts/ocl_test.hpp" + +#ifdef HAVE_OPENCL + +namespace cvtest { +namespace ocl { + +PARAM_TEST_CASE(FastNlMeansDenoisingTestBase, Channels, bool) +{ + int cn, templateWindowSize, searchWindowSize; + float h; + bool use_roi; + + TEST_DECLARE_INPUT_PARAMETER(src) + TEST_DECLARE_OUTPUT_PARAMETER(dst) + + virtual void SetUp() + { + cn = GET_PARAM(0); + use_roi = GET_PARAM(1); + + templateWindowSize = 7; + searchWindowSize = 21; + h = 3.0f; + } + + virtual void generateTestData() + { + const int type = CV_8UC(cn); + + Size roiSize = randomSize(1, MAX_VALUE); + Border srcBorder = randomBorder(0, use_roi ? MAX_VALUE : 0); + randomSubMat(src, src_roi, roiSize, srcBorder, type, 0, 255); + + Border dstBorder = randomBorder(0, use_roi ? MAX_VALUE : 0); + randomSubMat(dst, dst_roi, roiSize, dstBorder, type, 0, 255); + + UMAT_UPLOAD_INPUT_PARAMETER(src) + UMAT_UPLOAD_OUTPUT_PARAMETER(dst) + } +}; + +typedef FastNlMeansDenoisingTestBase FastNlMeansDenoising; + +OCL_TEST_P(FastNlMeansDenoising, Mat) +{ + for (int j = 0; j < test_loop_times; j++) + { + generateTestData(); + + OCL_OFF(cv::fastNlMeansDenoising(src_roi, dst_roi, h, templateWindowSize, searchWindowSize)); + OCL_ON(cv::fastNlMeansDenoising(usrc_roi, udst_roi, h, templateWindowSize, searchWindowSize)); + + OCL_EXPECT_MATS_NEAR(dst, 1) + } +} + +OCL_INSTANTIATE_TEST_CASE_P(Photo, FastNlMeansDenoising, Combine(Values((Channels)1), Bool())); + +} } // namespace cvtest::ocl + +#endif // HAVE_OPENCL