diff --git a/modules/ximgproc/doc/ximgproc.bib b/modules/ximgproc/doc/ximgproc.bib index 322b37121..a9b43c59b 100644 --- a/modules/ximgproc/doc/ximgproc.bib +++ b/modules/ximgproc/doc/ximgproc.bib @@ -77,3 +77,14 @@ year={2008}, organization={ACM} } + +@inproceedings{xu2011image, + title={Image smoothing via L 0 gradient minimization}, + author={Xu, Li and Lu, Cewu and Xu, Yi and Jia, Jiaya}, + booktitle={ACM Transactions on Graphics (TOG)}, + volume={30}, + number={6}, + pages={174}, + year={2011}, + organization={ACM} +} diff --git a/modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp b/modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp index 60be77765..f1030daf5 100644 --- a/modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp +++ b/modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp @@ -378,6 +378,19 @@ it should be 0.25. Setting it to 1.0 may lead to streaking artifacts. */ CV_EXPORTS_W void fastGlobalSmootherFilter(InputArray guide, InputArray src, OutputArray dst, double lambda, double sigma_color, double lambda_attenuation=0.25, int num_iter=3); +/** @brief Global image smoothing via L0 gradient minimization. + +@param src source image for filtering with unsigned 8-bit or signed 16-bit or floating-point depth. + +@param dst destination image. + +@param lambda parameter defining the smooth term weight. + +@param kappa parameter defining the increasing factor of the weight of the gradient data term. + +For more details about L0 Smoother, see the original paper @cite xu2011image. +*/ +CV_EXPORTS_W void l0Smooth(InputArray src, OutputArray dst, double lambda = 0.02, double kappa = 2.0); //! @} } } diff --git a/modules/ximgproc/perf/perf_l0_smooth.cpp b/modules/ximgproc/perf/perf_l0_smooth.cpp new file mode 100644 index 000000000..4ae7820f7 --- /dev/null +++ b/modules/ximgproc/perf/perf_l0_smooth.cpp @@ -0,0 +1,81 @@ +/* + * 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 L0SmoothTestParam; +typedef TestBaseWithParam L0SmoothTest; + +PERF_TEST_P(L0SmoothTest, perf, + Combine( + SZ_TYPICAL, + Values(CV_8U, CV_16U, CV_32F, CV_64F), + Values(1, 3)) +) +{ + L0SmoothTestParam params = GetParam(); + Size sz = get<0>(params); + int depth = get<1>(params); + int srcCn = get<2>(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()); + + RNG rnd(sz.height + depth + srcCn); + double lambda = rnd.uniform(0.01, 0.05); + double kappa = rnd.uniform(1.0, 3.0); + + TEST_CYCLE_N(1) + { + l0Smooth(src, dst, lambda, kappa); + } + + SANITY_CHECK_NOTHING(); +} +} diff --git a/modules/ximgproc/src/l0_smooth.cpp b/modules/ximgproc/src/l0_smooth.cpp new file mode 100644 index 000000000..c4c30cda3 --- /dev/null +++ b/modules/ximgproc/src/l0_smooth.cpp @@ -0,0 +1,308 @@ +/* + * 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 "precomp.hpp" +#include +#include +#include + +using namespace cv; +using namespace std; + +namespace +{ + void shift(InputArray src, OutputArray dst, int shift_x, int shift_y) { + Mat S = src.getMat(); + Mat D = dst.getMat(); + + if(S.data == D.data){ + S = S.clone(); + } + + D.create(S.size(), S.type()); + + Mat s0(S, Rect(0, 0, S.cols - shift_x, S.rows - shift_y)); + Mat s1(S, Rect(S.cols - shift_x, 0, shift_x, S.rows - shift_y)); + Mat s2(S, Rect(0, S.rows - shift_y, S.cols-shift_x, shift_y)); + Mat s3(S, Rect(S.cols - shift_x, S.rows- shift_y, shift_x, shift_y)); + + Mat d0(D, Rect(shift_x, shift_y, S.cols - shift_x, S.rows - shift_y)); + Mat d1(D, Rect(0, shift_y, shift_x, S.rows - shift_y)); + Mat d2(D, Rect(shift_x, 0, S.cols-shift_x, shift_y)); + Mat d3(D, Rect(0,0,shift_x, shift_y)); + + s0.copyTo(d0); + s1.copyTo(d1); + s2.copyTo(d2); + s3.copyTo(d3); + } + + // dft after padding imaginary + void fft(InputArray src, OutputArray dst) { + Mat S = src.getMat(); + Mat planes[] = { S, Mat::zeros(S.size(), S.type()) }; + merge(planes, 2, dst); + + // compute the result + dft(dst, dst); + } + + void psf2otf(InputArray src, OutputArray dst, int height, int width){ + Mat S = src.getMat(); + Mat D = dst.getMat(); + + Mat padded; + + if(S.data == D.data){ + S = S.clone(); + } + + // add padding + copyMakeBorder(S, padded, 0, height - S.rows, 0, width - S.cols, + BORDER_CONSTANT, Scalar::all(0)); + + shift(padded, padded, width - S.cols / 2, height - S.rows / 2); + + // convert to frequency domain + fft(padded, dst); + } + + void dftMultiChannel(InputArray src, vector &dst){ + Mat S = src.getMat(); + + split(S, dst); + + for(int i = 0; i < S.channels(); i++){ + fft(dst[i], dst[i]); + } + } + + void idftMultiChannel(const vector &src, OutputArray dst){ + Mat *channels = new Mat[src.size()]; + + for(int i = 0 ; unsigned(i) < src.size(); i++){ + idft(src[i], channels[i]); + Mat realImg[2]; + split(channels[i], realImg); + channels[i] = realImg[0] / src[i].cols / src[i].rows; + } + + Mat D; + merge(channels, src.size(), D); + D.copyTo(dst); + + delete[] channels; + } + + void addComplex(InputArray aSrc, int bSrc, OutputArray dst){ + Mat panels[2]; + split(aSrc.getMat(), panels); + panels[0] = panels[0] + bSrc; + merge(panels, 2, dst); + } + + void divComplexByReal(InputArray aSrc, InputArray bSrc, OutputArray dst){ + Mat aPanels[2]; + Mat bPanels[2]; + split(aSrc.getMat(), aPanels); + split(bSrc.getMat(), bPanels); + + Mat realPart; + Mat imaginaryPart; + + divide(aPanels[0], bSrc.getMat(), realPart); + divide(aPanels[1], bSrc.getMat(), imaginaryPart); + + aPanels[0] = realPart; + aPanels[1] = imaginaryPart; + + Mat rst; + merge(aPanels, 2, dst); + } + + void divComplexByRealMultiChannel(const vector &numer, + const vector &denom, vector &dst) + { + for(int i = 0; unsigned(i) < numer.size(); i++) + { + divComplexByReal(numer[i], denom[i], dst[i]); + } + } + + // power of 2 of the absolute value of the complex + Mat pow2absComplex(InputArray src){ + Mat S = src.getMat(); + + Mat sPanels[2]; + split(S, sPanels); + + return sPanels[0].mul(sPanels[0]) + sPanels[1].mul(sPanels[1]); + } +} + +namespace cv +{ +namespace ximgproc +{ + +void l0Smooth(InputArray src, OutputArray dst, double lambda, double kappa) +{ + Mat S = src.getMat(); + + CV_Assert(!S.empty()); + CV_Assert(S.depth() == CV_8U || S.depth() == CV_16U + || S.depth() == CV_32F || S.depth() == CV_64F); + + dst.create(src.size(), src.type()); + + if(S.data == dst.getMat().data){ + S = S.clone(); + } + + if(S.depth() == CV_8U) + { + S.convertTo(S, CV_32F, 1/255.0f); + } + else if(S.depth() == CV_16U) + { + S.convertTo(S, CV_32F, 1/65535.0f); + }else if(S.depth() == CV_64F){ + S.convertTo(S, CV_32F); + } + + const double betaMax = 100000; + + // gradient operators in frequency domain + Mat otfFx, otfFy; + float kernel[2] = {-1, 1}; + float kernel_inv[2] = {1,-1}; + psf2otf(Mat(1,2,CV_32FC1, kernel_inv), otfFx, S.rows, S.cols); + psf2otf(Mat(2,1,CV_32FC1, kernel_inv), otfFy, S.rows, S.cols); + + vector denomConst; + Mat tmp = pow2absComplex(otfFx) + pow2absComplex(otfFy); + + for(int i = 0; i < S.channels(); i++){ + denomConst.push_back(tmp); + } + + // input image in frequency domain + vector numerConst; + dftMultiChannel(S, numerConst); + + /********************************* + * solver + *********************************/ + double beta = 2 * lambda; + while(beta < betaMax){ + // h, v subproblem + Mat h, v; + + filter2D(S, h, -1, Mat(1, 2, CV_32FC1, kernel), Point(0, 0), + 0, BORDER_REPLICATE); + filter2D(S, v, -1, Mat(2, 1, CV_32FC1, kernel), Point(0, 0), + 0, BORDER_REPLICATE); + + Mat hvMag = h.mul(h) + v.mul(v); + + Mat mask; + if(S.channels() == 1) + { + threshold(hvMag, mask, lambda/beta, 1, THRESH_BINARY); + } + else if(S.channels() > 1) + { + Mat *channels = new Mat[S.channels()]; + split(hvMag, channels); + hvMag = channels[0]; + + for(int i = 1; i < S.channels(); i++){ + hvMag = hvMag + channels[i]; + } + + threshold(hvMag, mask, lambda/beta, 1, THRESH_BINARY); + + Mat in[] = {mask, mask, mask}; + merge(in, 3, mask); + + delete[] channels; + } + + h = h.mul(mask); + v = v.mul(mask); + + // S subproblem + vector denom(S.channels()); + for(int i = 0; i < S.channels(); i++){ + denom[i] = beta * denomConst[i] + 1; + } + + Mat hGrad, vGrad; + filter2D(h, hGrad, -1, Mat(1, 2, CV_32FC1, kernel_inv)); + filter2D(v, vGrad, -1, Mat(2, 1, CV_32FC1, kernel_inv)); + + vector hvGradFreq; + dftMultiChannel(hGrad+vGrad, hvGradFreq); + + vector numer(S.channels()); + for(int i = 0; i < S.channels(); i++){ + numer[i] = numerConst[i] + hvGradFreq[i] * beta; + } + + vector sFreq(S.channels()); + divComplexByRealMultiChannel(numer, denom, sFreq); + + idftMultiChannel(sFreq, S); + + beta = beta * kappa; + } + + Mat D = dst.getMat(); + if(D.depth() == CV_8U) + { + S.convertTo(D, CV_8U, 255); + } + else if(D.depth() == CV_16U) + { + S.convertTo(D, CV_16U, 65535); + }else if(D.depth() == CV_64F){ + S.convertTo(D, CV_64F); + }else{ + S.copyTo(D); + } +} +} +} diff --git a/modules/ximgproc/test/test_l0_smooth.cpp b/modules/ximgproc/test/test_l0_smooth.cpp new file mode 100644 index 000000000..dae1b4341 --- /dev/null +++ b/modules/ximgproc/test/test_l0_smooth.cpp @@ -0,0 +1,120 @@ +/* + * 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; + +CV_ENUM(SrcTypes, CV_8UC1, CV_8UC3, CV_16UC1, CV_16UC3); +typedef tuple L0SmoothParams; +typedef TestWithParam L0SmoothTest; + +TEST(L0SmoothTest, SplatSurfaceAccuracy) +{ + RNG rnd(0); + + for (int i = 0; i < 3; i++) + { + Size sz(rnd.uniform(512, 1024), rnd.uniform(512, 1024)); + + Scalar surfaceValue; + int srcCn = 3; + rnd.fill(surfaceValue, RNG::UNIFORM, 0, 255); + Mat src(sz, CV_MAKE_TYPE(CV_8U, srcCn), surfaceValue); + + double lambda = rnd.uniform(0.01, 0.05); + double kappa = rnd.uniform(1.5, 5.0); + + Mat res; + l0Smooth(src, res, lambda, kappa); + + // When filtering a constant image we should get the same image: + double normL1 = cvtest::norm(src, res, NORM_L1)/src.total()/src.channels(); + EXPECT_LE(normL1, 1.0/64); + } +} + +TEST_P(L0SmoothTest, MultiThreadReproducibility) +{ + if (cv::getNumberOfCPUs() == 1) + return; + + double MAX_DIF = 10.0; + double MAX_MEAN_DIF = 1.0 / 8.0; + int loopsCount = 2; + RNG rng(0); + + L0SmoothParams params = GetParam(); + Size size = get<0>(params); + int srcType = get<1>(params); + + Mat src(size,srcType); + if(src.depth()==CV_8U) + randu(src, 0, 255); + else if(src.depth()==CV_16U) + randu(src, 0, 65535); + else + randu(src, -100000.0f, 100000.0f); + + + for (int iter = 0; iter <= loopsCount; iter++) + { + double lambda = rng.uniform(0.01, 0.05); + double kappa = rng.uniform(1.5, 5.0); + + cv::setNumThreads(cv::getNumberOfCPUs()); + Mat resMultiThread; + l0Smooth(src, resMultiThread, lambda, kappa); + + cv::setNumThreads(1); + Mat resSingleThread; + l0Smooth(src, resSingleThread, lambda, kappa); + + 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(FullSet, L0SmoothTest,Combine(Values(szODD, szQVGA), SrcTypes::all())); + +}