Merge pull request #382 from zhou-chao:l0smooth

modules/ximgproc/doc/ximgproc.bib

@@ -89,6 +89,18 @@
   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}
+}
+
 @inproceedings{Revaud2015,
   title={EpicFlow: Edge-Preserving Interpolation of Correspondences for Optical Flow},
   author={Revaud, Jerome and Weinzaepfel, Philippe and Harchaoui, Zaid and Schmid, Cordelia},

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);
 //! @}
 }
 }

modules/ximgproc/perf/perf_l0_smooth.cpp

+/*
+ *  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<Size, MatType, int> L0SmoothTestParam;
+typedef TestBaseWithParam<L0SmoothTestParam> 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();
+}
+}

modules/ximgproc/src/l0_smooth.cpp

+/*
+ *  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 <vector>
+#include <opencv2/core.hpp>
+#include <opencv2/imgproc.hpp>
+
+using namespace cv;
+using namespace std;
+
+namespace
+{
+
+    class ParallelDft : public ParallelLoopBody
+    {
+    private:
+        vector<Mat> src_;
+    public:
+        ParallelDft(vector<Mat> &s)
+        {
+            src_ = s;
+        }
+        void operator() (const Range& range) const
+        {
+            for (int i = range.start; i != range.end; i++)
+            {
+                dft(src_[i], src_[i]);
+            }
+        }
+    };
+
+    class ParallelIdft : public ParallelLoopBody
+    {
+    private:
+        vector<Mat> src_;
+    public:
+        ParallelIdft(vector<Mat> &s)
+        {
+            src_ = s;
+        }
+        void operator() (const Range& range) const
+        {
+            for (int i = range.start; i != range.end; i++)
+            {
+                idft(src_[i], src_[i],DFT_SCALE);
+            }
+        }
+    };
+
+    class ParallelDivComplexByReal : public ParallelLoopBody
+    {
+    private:
+        vector<Mat> numer_;
+        vector<Mat> denom_;
+        vector<Mat> dst_;
+
+    public:
+        ParallelDivComplexByReal(vector<Mat> &numer, vector<Mat> &denom, vector<Mat> &dst)
+        {
+            numer_ = numer;
+            denom_ = denom;
+            dst_ = dst;
+        }
+        void operator() (const Range& range) const
+        {
+            for (int i = range.start; i != range.end; i++)
+            {
+                Mat aPanels[2];
+                Mat bPanels[2];
+                split(numer_[i], aPanels);
+                split(denom_[i], bPanels);
+
+                Mat realPart;
+                Mat imaginaryPart;
+
+                divide(aPanels[0], denom_[i], realPart);
+                divide(aPanels[1], denom_[i], imaginaryPart);
+
+                aPanels[0] = realPart;
+                aPanels[1] = imaginaryPart;
+
+                merge(aPanels, 2, dst_[i]);
+            }
+        }
+    };
+
+
+    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.clone(), Mat::zeros(S.size(), S.type())};
+        Mat x;
+        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<Mat> &dst)
+    {
+        Mat S = src.getMat();
+
+        split(S, dst);
+
+        for(int i = 0; i < S.channels(); i++){
+            Mat planes[] = {dst[i].clone(), Mat::zeros(dst[i].size(), dst[i].type())};
+            merge(planes, 2, dst[i]);
+        }
+
+        parallel_for_(cv::Range(0,S.channels()), ParallelDft(dst));
+
+    }
+
+    void idftMultiChannel(const vector<Mat> &src, OutputArray dst)
+    {
+        vector<Mat> channels(src);
+
+        parallel_for_(Range(0, int(src.size())), ParallelIdft(channels));
+
+        for(int i = 0; unsigned(i) < src.size(); i++){
+            Mat panels[2];
+            split(channels[i], panels);
+            channels[i] = panels[0];
+        }
+
+        Mat D;
+        merge(channels, D);
+        D.copyTo(dst);
+    }
+
+    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 divComplexByRealMultiChannel(vector<Mat> &numer,
+        vector<Mat> &denom, vector<Mat> &dst)
+    {
+
+        for(int i = 0; unsigned(i) < numer.size(); i++)
+        {
+            dst[i].create(numer[i].size(), numer[i].type());
+        }
+        parallel_for_(Range(0, int(numer.size())), ParallelDivComplexByReal(numer, denom, dst));
+
+    }
+
+        // power of 2 of the absolute value of the complex
+    Mat pow2absComplex(InputArray src)
+    {
+        Mat S = src.getMat();
+
+        Mat sPanels[2];
+        split(S, sPanels);
+
+        Mat mag;
+        magnitude(sPanels[0], sPanels[1], mag);
+        pow(mag, 2, mag);
+
+        return mag;
+    }
+}
+
+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<Mat> 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<Mat> 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)
+                {
+                    vector<Mat> channels(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);
+                }
+
+                h = h.mul(mask);
+                v = v.mul(mask);
+
+                // S subproblem
+                vector<Mat> 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<Mat> hvGradFreq;
+                dftMultiChannel(hGrad+vGrad, hvGradFreq);
+
+                vector<Mat> numer(S.channels());
+                for(int i = 0; i < S.channels(); i++)
+                {
+                    numer[i] = numerConst[i] + hvGradFreq[i] * beta;
+                }
+
+                vector<Mat> 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);
+            }
+        }
+    }
+}

modules/ximgproc/test/test_l0_smooth.cpp

+/*
+ *  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<Size, SrcTypes> L0SmoothParams;
+typedef TestWithParam<L0SmoothParams> 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()));
+
+}