diff --git a/modules/ximgproc/doc/ximgproc.bib b/modules/ximgproc/doc/ximgproc.bib index 0fbbb374b..06060b60c 100644 --- a/modules/ximgproc/doc/ximgproc.bib +++ b/modules/ximgproc/doc/ximgproc.bib @@ -149,6 +149,17 @@ year = {2015} } +@article{Cho2014, + title = {Bilateral Texture Filtering}, + author = {Hojin Cho and Hyunjoon Lee and Henry Kang and Seungyong Lee}, + journal = {ACM Transactions on Graphics}, + year = {2014}, + volume = {33}, + number = {4}, + month = {July}, + pages = {128:1--128:8} +} + @incollection{zhang2014rolling, title={Rolling guidance filter}, author={Zhang, Qi and Shen, Xiaoyong and Xu, Li and Jia, Jiaya}, diff --git a/modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp b/modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp index 9b722fa75..4709798b6 100644 --- a/modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp +++ b/modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp @@ -316,6 +316,28 @@ proportional to sigmaSpace . CV_EXPORTS_W void jointBilateralFilter(InputArray joint, InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace, int borderType = BORDER_DEFAULT); +/** @brief Applies the bilateral texture filter to an image. It performs structure-preserving texture filter. +For more details about this filter see @cite Cho2014. + +@param src Source image whose depth is 8-bit UINT or 32-bit FLOAT + +@param dst Destination image of the same size and type as src. + +@param fr Radius of kernel to be used for filtering. It should be positive integer + +@param numIter Number of iterations of algorithm, It should be positive integer + +@param sigmaAlpha Controls the sharpness of the weight transition from edges to smooth/texture regions, where +a bigger value means sharper transition. When the value is negative, it is automatically calculated. + +@param sigmaAvg Range blur parameter for texture blurring. Larger value makes result to be more blurred. When the +value is negative, it is automatically calculated as described in the paper. + +@sa rollingGuidanceFilter, bilateralFilter +*/ +CV_EXPORTS_W +void bilateralTextureFilter(InputArray src, OutputArray dst, int fr = 3, int numIter = 1, double sigmaAlpha = -1., double sigmaAvg = -1.); + ////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////// diff --git a/modules/ximgproc/perf/perf_bilateral_texture_filter.cpp b/modules/ximgproc/perf/perf_bilateral_texture_filter.cpp new file mode 100644 index 000000000..0da51babe --- /dev/null +++ b/modules/ximgproc/perf/perf_bilateral_texture_filter.cpp @@ -0,0 +1,83 @@ +/* + * 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 + * (3 - clause BSD License) + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met : + * + * *Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * * Redistributions 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. + * + * * Neither the names of the copyright holders nor the names of the contributors + * may 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 copyright holders 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. + */ + +#include "perf_precomp.hpp" + +namespace cvtest +{ + +using std::tr1::tuple; +using std::tr1::get; +using namespace perf; +using namespace testing; +using namespace cv; +using namespace cv::ximgproc; + +typedef tuple BTFTestParam; +typedef TestBaseWithParam BilateralTextureFilterTest; + +PERF_TEST_P(BilateralTextureFilterTest, perf, + Combine( + Values(2), + Values(0.5), + Values(0.5), + SZ_TYPICAL, + Values(CV_8U, CV_32F), + Values(1, 3)) +) +{ + BTFTestParam params = GetParam(); + int fr = get<0>(params); + double sigmaAlpha = get<1>(params); + double sigmaAvg = get<2>(params); + Size sz = get<3>(params); + int depth = get<4>(params); + int srcCn = get<5>(params); + + Mat src(sz, CV_MAKE_TYPE(depth,srcCn)); + Mat dst(sz, src.type()); + + cv::setNumThreads(cv::getNumberOfCPUs()); + declare.in(src, WARMUP_RNG).out(dst).tbb_threads(cv::getNumberOfCPUs()); + + TEST_CYCLE_N(1) + { + bilateralTextureFilter(src, dst, fr, 1, sigmaAlpha, sigmaAvg); + } + + SANITY_CHECK_NOTHING(); +} +} diff --git a/modules/ximgproc/src/bilateral_texture_filter.cpp b/modules/ximgproc/src/bilateral_texture_filter.cpp new file mode 100644 index 000000000..5be6c0a6b --- /dev/null +++ b/modules/ximgproc/src/bilateral_texture_filter.cpp @@ -0,0 +1,357 @@ +/*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-2011, Willow Garage Inc., all rights reserved. +// Third party copyrights are property of their respective owners. +// +// * 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 Intel Corporation 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 +#include + +namespace cv +{ +namespace ximgproc +{ + void compute_mRTV(const Mat& L, Mat& mRTV, int fr); + void compute_G(const Mat& B, const Mat& mRTV, Mat& G, Mat& alpha, int fr); + void joint_bilateral_filter(const Mat& img, const Mat& G, Mat& r_img, int fr2, double sigma_avg); + void joint_bilateral_filter3(const Mat& img, const Mat& G, Mat& r_img, int fr2, double sigma_avg); + + void bilateralTextureFilter(InputArray src_, OutputArray dst_, int fr, + int numIter, double sigmaAlpha, double sigmaAvg) + { + CV_Assert(!src_.empty()); + + Mat src = src_.getMat(); + CV_Assert(src.depth() == CV_8U || src.depth() == CV_32F); + + CV_Assert(fr > 0 && numIter > 0); + + if (sigmaAlpha < 0) + sigmaAlpha = 5. * fr; + if (sigmaAvg < 0) + sigmaAvg = 0.05 * sqrt(src.channels()); + + Mat I; + src.copyTo(I); + if (src.type() == CV_8UC1) { + I.convertTo(I, CV_32FC1, 1.0 / 255.0); + } + else if (src.type() == CV_8UC3) { + I.convertTo(I, CV_32FC3, 1.0 / 255.0); + } + + for (int iter = 0; iter < numIter; iter++) + { + Mat B; + blur(I, B, Size(2 * fr + 1, 2 * fr + 1), Point(-1, -1), BORDER_REFLECT); + + Mat mRTV; + compute_mRTV(I, mRTV, fr); + + Mat G, minmRTV; + compute_G(B, mRTV, G, minmRTV, fr); + + // alpha blending + Mat Gtilde; + Mat diff = mRTV - minmRTV; + Mat alpha = -diff.mul(sigmaAlpha); + exp(alpha, alpha); + alpha = alpha + 1.; + pow(alpha, -1, alpha); + alpha = (alpha - 0.5) * 2; + Mat alphainv = -(alpha - 1); + + std::vector Gi, Bi; + Gi.resize(I.channels()); + Bi.resize(I.channels()); + if (I.channels() == 3) { + split(G, &Gi[0]); + split(B, &Bi[0]); + } + else { + G.copyTo(Gi[0]); + B.copyTo(Bi[0]); + } + + std::vector Gtildei; + Gtildei.resize(I.channels()); + for (int i = 0; i < B.channels(); i++) + Gtildei[i] = Gi[i].mul(alpha) + Bi[i].mul(alphainv); + merge(&Gtildei[0], B.channels(), Gtilde); + + // joint bilateral filter + cv::Mat J; + if (I.channels() == 1) + joint_bilateral_filter(I, Gtilde, J, fr * 2, sigmaAvg); + else if (I.channels() == 3) + joint_bilateral_filter3(I, Gtilde, J, fr * 2, sigmaAvg); + I = J; + } + if (src.type() == CV_8UC1) { + I.convertTo(I, CV_8UC1, 255.0); + } + else if (src.type() == CV_8UC3) { + I.convertTo(I, CV_8UC3, 255.0); + } + + I.copyTo(dst_); + } + + void compute_mRTV(const Mat& L, Mat& mRTV, int fr) + { + mRTV = Mat::zeros(L.size(), CV_32FC1); + + const float eps = 0.00001f; + + // Calculate image derivative(gradient) + Mat G; + Mat Gx, Gy, kernelx, kernely; + kernelx = Mat::zeros(1, 3, CV_32F); + kernelx.at(0, 1) = -1.0; + kernelx.at(0, 2) = 1.0; + filter2D(L, Gx, -1, kernelx, Point(-1, -1), 0, BORDER_REFLECT); + kernely = Mat::zeros(3, 1, CV_32F); + kernely.at(1, 0) = -1.0; + kernely.at(2, 0) = 1.0; + filter2D(L, Gy, -1, kernely, Point(-1, -1), 0, BORDER_REFLECT); + + Gx = Gx.mul(Gx); + Gy = Gy.mul(Gy); + sqrt(Gx + Gy, G); + + // Pad image L and G + Mat padL; + Mat padG; + copyMakeBorder(L, padL, fr, fr, fr, fr, BORDER_REFLECT); + copyMakeBorder(G, padG, fr, fr, fr, fr, BORDER_REFLECT); + + // Calculate maxL, minL, maxG, sumG + int pu = fr; + int pb = pu + L.rows; + int pl = fr; + int pr = pl + L.cols; + + std::vector Li, Gi; + Li.resize(L.channels()); + Gi.resize(L.channels()); + if (L.channels() == 3) { + split(padL, &Li[0]); + split(padG, &Gi[0]); + } + else { + padL.copyTo(Li[0]); + padG.copyTo(Gi[0]); + } + + for (int i = 0; i < L.channels(); i++) + { + Mat maxL = Mat::zeros(L.size(), CV_32FC1); + Mat minL = Mat::ones(L.size(), CV_32FC1); + Mat maxG = Mat::zeros(L.size(), CV_32FC1); + Mat sumG = Mat::zeros(L.size(), CV_32FC1); + for (int y = -fr; y <= fr; y++) + { + for (int x = -fr; x <= fr; x++) + { + Mat temp = Li[i]( + Range(pu + y, pb + y), + Range(pl + x, pr + x) + ); + maxL = max(maxL, temp); + minL = min(minL, temp); + + temp = Gi[i]( + Range(pu + y, pb + y), + Range(pl + x, pr + x) + ); + maxG = max(maxG, temp); + sumG = sumG + temp; + } + } + Mat deltai = maxL - minL; + sumG = max(sumG, eps); + Mat mRTVi = maxG / sumG * (2 * fr + 1); + mRTV = mRTV + mRTVi.mul(deltai); + } + if (L.channels() == 3) + mRTV = mRTV / 3; + } + + void compute_G(const Mat& B, const Mat& mRTV, Mat& G, Mat& alpha, int fr) + { + B.copyTo(G); + alpha = Mat::ones(B.size(), CV_32FC1); + for (int y = -fr; y <= fr; y++) + { + for (int x = -fr; x <= fr; x++) + { + Point pb; + Point pt; + for (pb.y = 0; pb.y < B.rows; pb.y++) + { + for (pb.x = 0; pb.x < B.cols; pb.x++) + { + pt.x = min(max(pb.x + x, 0), B.cols - 1); + pt.y = min(max(pb.y + y, 0), B.rows - 1); + if (alpha.at(pb) > mRTV.at(pt)) + { + alpha.at(pb) = mRTV.at(pt); + if (B.channels() == 3) + G.at(pb) = B.at(pt); + else if (B.channels() == 1) + G.at(pb) = B.at(pt); + } + } + } + } + } + } + + void joint_bilateral_filter(const Mat& img, const Mat& G, Mat& r_img, int fr2, double sigma_avg) + { + Mat p_G; + copyMakeBorder(G, p_G, fr2, fr2, fr2, fr2, BORDER_REFLECT); + + Mat p_img; + copyMakeBorder(img, p_img, fr2, fr2, fr2, fr2, BORDER_REFLECT); + + Mat SW; + if (SW.empty()) { + SW = Mat(2*fr2+1, 2*fr2+1, CV_32FC1); + int r, c; + float y, x; + for (r = 0, y = (float)-fr2; r < SW.rows; r++, y += 1.0) { + for(c = 0, x = (float)-fr2; c < SW.cols; c++, x += 1.0) { + SW.at(r,c) = exp(-(x*x + y*y) / (2*fr2*fr2)); + } + } + } + + r_img = Mat::zeros(img.size(), CV_32FC1); + { + Mat sum_d_W = Mat::zeros(img.size(), CV_32FC1); + Mat d_W = Mat::zeros(G.size(), CV_32FC1); + + for (int x = -fr2; x <= fr2; x++) { + for (int y = -fr2; y <= fr2; y++) { + d_W = p_G(Rect(fr2+x, fr2+y, img.cols, img.rows)) - G; + multiply(d_W, d_W, d_W); + exp(-0.5 * d_W / (sigma_avg*sigma_avg), d_W); + + d_W = d_W * SW.at(fr2+y, fr2+x); //Gaussian weight + + sum_d_W = sum_d_W + d_W; + multiply(d_W, p_img(Rect(fr2+x, fr2+y, img.cols, img.rows)), d_W); + r_img = r_img + d_W; + } + } + max(1e-5f, sum_d_W, sum_d_W); + divide(r_img, sum_d_W, r_img); + } + } + + void joint_bilateral_filter3(const Mat& img, const Mat& G, Mat& r_img, int fr2, double sigma_avg) + { + Mat p_G; + copyMakeBorder(G, p_G, fr2, fr2, fr2, fr2, BORDER_REFLECT); + + Mat p_img; + copyMakeBorder(img, p_img, fr2, fr2, fr2, fr2, BORDER_REFLECT); + + Mat SW; + if (SW.empty()) { + SW = Mat(2*fr2+1, 2*fr2+1, CV_32FC1); + int r, c; + float y, x; + for (r = 0, y = (float)-fr2; r < SW.rows; r++, y += 1.0) { + for(c = 0, x = (float)-fr2; c < SW.cols; c++, x += 1.0) { + SW.at(r,c) = exp(-(x*x + y*y) / (2*fr2*fr2)); + } + } + } + + std::vector G_channels(3); + split(G, G_channels); + std::vector p_G_channels(3); + split(p_G, p_G_channels); + std::vector p_img_channels(3); + split(p_img, p_img_channels); + + Mat sum_d_W = Mat::zeros(img.size(), CV_32FC1); + std::vector d_W_channels(3); + for (int ch = 0; ch < 3; ch++) { + d_W_channels[ch] = Mat::zeros(G.size(), CV_32FC1); + } + Mat d_W = Mat::zeros(G.size(), CV_32FC1); + + std::vector r_img_channels(3); + for (int ch = 0; ch < 3; ch++) { + r_img_channels[ch] = Mat::zeros(G.size(), CV_32FC1); + } + + for (int x = -fr2; x <= fr2; x++) { + for (int y = -fr2; y <= fr2; y++) { + d_W.setTo(0); + for (int ch = 0; ch < 3; ch++) { + subtract(p_G_channels[ch](Rect(fr2+x, fr2+y, img.cols, img.rows)), G_channels[ch], d_W_channels[ch]); + multiply(d_W_channels[ch], d_W_channels[ch], d_W_channels[ch]); + } + for (int ch = 0; ch < 3; ch++) { + add(d_W, d_W_channels[ch], d_W); + } + exp(-0.5 * d_W / (sigma_avg*sigma_avg), d_W); + + d_W = d_W * SW.at(fr2+y, fr2+x); //Gaussian weight + + add(sum_d_W, d_W, sum_d_W); + + for (int ch = 0; ch < 3; ch++) { + Mat n_p_img = p_img_channels[ch](Rect(fr2+x, fr2+y, img.cols, img.rows)); + accumulateProduct(d_W, n_p_img, r_img_channels[ch]); + } + } + } + + max(1e-5f, sum_d_W, sum_d_W); + for (int ch = 0; ch < 3; ch++) { + divide(r_img_channels[ch], sum_d_W, r_img_channels[ch]); + } + merge(r_img_channels, r_img); + } + +} +} diff --git a/modules/ximgproc/test/test_bilateral_texture_filter.cpp b/modules/ximgproc/test/test_bilateral_texture_filter.cpp new file mode 100644 index 000000000..2c53e2377 --- /dev/null +++ b/modules/ximgproc/test/test_bilateral_texture_filter.cpp @@ -0,0 +1,141 @@ +/* + * 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 + * (3 - clause BSD License) + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met : + * + * *Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * * Redistributions 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. + * + * * Neither the names of the copyright holders nor the names of the contributors + * may 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 copyright holders 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. + */ + +#include "test_precomp.hpp" + +namespace cvtest +{ + +using namespace std; +using namespace std::tr1; +using namespace testing; +using namespace perf; +using namespace cv; +using namespace cv::ximgproc; + +////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////// + +typedef tuple BTFParams; +typedef TestWithParam BilateralTextureFilterTest; + +TEST_P(BilateralTextureFilterTest, SplatSurfaceAccuracy) +{ + BTFParams params = GetParam(); + int fr = get<0>(params); + double sigmaAlpha = get<1>(params); + double sigmaAvg = get<2>(params); + int depth = get<3>(params); + int srcCn = get<4>(params); + + RNG rnd(0); + + Size sz(rnd.uniform(256,512), rnd.uniform(256,512)); + + for (int i = 0; i < 5; i++) + { + Scalar surfaceValue; + if(depth == CV_8U) + rnd.fill(surfaceValue, RNG::UNIFORM, 0, 255); + else + rnd.fill(surfaceValue, RNG::UNIFORM, 0.0f, 1.0f); + + Mat src(sz, CV_MAKE_TYPE(depth, srcCn), surfaceValue); + + Mat res; + bilateralTextureFilter(src, res, fr, 1, sigmaAlpha, sigmaAvg); + + double normL1 = cvtest::norm(src, res, NORM_L1)/src.total()/src.channels(); + EXPECT_LE(normL1, 1.0/64.0); + } +} + +TEST_P(BilateralTextureFilterTest, MultiThreadReproducibility) +{ + if (cv::getNumberOfCPUs() == 1) + return; + + BTFParams params = GetParam(); + int fr = get<0>(params); + double sigmaAlpha = get<1>(params); + double sigmaAvg = get<2>(params); + int depth = get<3>(params); + int srcCn = get<4>(params); + + double MAX_DIF = 1.0; + double MAX_MEAN_DIF = 1.0 / 64.0; + int loopsCount = 2; + RNG rnd(1); + + Size sz(rnd.uniform(256,512), rnd.uniform(256,512)); + + Mat src(sz,CV_MAKE_TYPE(depth, srcCn)); + if(src.depth()==CV_8U) + randu(src, 0, 255); + else if(src.depth()==CV_16S) + randu(src, -32767, 32767); + else + randu(src, 0.0f, 1.0f); + + for (int iter = 0; iter <= loopsCount; iter++) + { + cv::setNumThreads(cv::getNumberOfCPUs()); + Mat resMultiThread; + bilateralTextureFilter(src, resMultiThread, fr, 1, sigmaAlpha, sigmaAvg); + + cv::setNumThreads(1); + Mat resSingleThread; + bilateralTextureFilter(src, resSingleThread, fr, 1, sigmaAlpha, sigmaAvg); + + EXPECT_LE(cv::norm(resSingleThread, resMultiThread, NORM_INF), MAX_DIF); + EXPECT_LE( + cv::norm(resSingleThread, resMultiThread, NORM_L1), + MAX_MEAN_DIF * src.total() * src.channels()); + } +} + +INSTANTIATE_TEST_CASE_P( + TypicalSet1, + BilateralTextureFilterTest, + Combine( + Values(2), + Values(0.5), + Values(0.5), + Values(CV_8U, CV_32F), + Values(1, 3) + ) +); +}