Merge pull request #1819 from berak:fast_bilateral

The Fast Bilateral Filter (#1819)

PR is created by: berak <px1704@web.de>
based on work https://github.com/opencv/opencv_contrib/pull/1317

* fbs_filter v1.0 has been contributed

* use boost unordered_map

* add brief description for fbs_filter

* fix format

* fix channels bug

* modify doc for fbs_filter

* check c++ 11

* asDiagonal -> diagonal

* rosolve warning

* fix eigen3 dependency

* Eigen/Core

* test HEAV_EIGEN

* setZero bug

* unordered_map test

* fix macro bug

* fix boost not found

* fix eigen macro bug

* fix eigen macro bug

* fix eigen macro bug

* fix eigen macro bug

*  add test file

*  fix test macro

*  fix test macro

*   add test

*   add test

* add sample colorize

* fix macro

*  fix colorize.cpp

*  fix colorize.cpp

*  fix colorize.cpp

*  fix colorize.cpp

* add fbs filter demo

* add fbs filter demo

* add fbs filter demo

* use fgsfilter for guess

* add parameter num_iter and max_tol

* add a option for colorize sample

* add  bibtex

* add  bibtex

* fix a colorize demo bug

* size optimize

* taking over the fast bilateral solver

* taking over the fast bilateral solver

* try to fix test_fbs_filter

* missed a bib bracket

modules/ximgproc/doc/ximgproc.bib

@@ -309,5 +309,11 @@
   title = {Ridge Filter Mathematica},
   author = {Wolfram Mathematica},
   url = {http://reference.wolfram.com/language/ref/RidgeFilter.html}
+}
 
+@article{BarronPoole2016,
+  author = {Jonathan T Barron and Ben Poole},
+  title = {The Fast Bilateral Solver},
+  journal = {ECCV},
+  year = {2016},
 }

modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp

@@ -375,6 +375,67 @@ void rollingGuidanceFilter(InputArray src, OutputArray dst, int d = -1, double s
 //////////////////////////////////////////////////////////////////////////
 //////////////////////////////////////////////////////////////////////////
 
+/** @brief Interface for implementations of The Fast Bilateral Solver.
+
+For more details about this solver see @cite BarronPoole2016 .
+*/
+class CV_EXPORTS_W FastBilateralSolverFilter : public Algorithm
+{
+public:
+    /** @brief Apply smoothing operation to the source image.
+
+    @param src source image for filtering with unsigned 8-bit or signed 16-bit or floating-point 32-bit depth and up to 4 channels.
+
+    @param confidence confidence image with unsigned 8-bit or signed 16-bit or floating-point 32-bit confidence and 1 channel.
+
+    @param dst destination image.
+    */
+    CV_WRAP virtual void filter(InputArray src, InputArray confidence, OutputArray dst) = 0;
+};
+
+/** @brief Factory method, create instance of FastBilateralSolverFilter and execute the initialization routines.
+
+@param guide image serving as guide for filtering. It should have 8-bit depth and either 1 or 3 channels.
+
+@param sigma_spatial parameter, that is similar to spatial space sigma in bilateralFilter.
+
+@param sigma_luma parameter, that is similar to luma space sigma in bilateralFilter.
+
+@param sigma_chroma parameter, that is similar to chroma space sigma in bilateralFilter.
+
+@param num_iter number of iterations used for solving, 25 is usually enough.
+
+@param max_tol solving tolerance used for solving, 25 is usually enough.
+
+For more details about the Fast Bilateral Solver parameters, see the original paper @cite BarronPoole2016.
+*/
+CV_EXPORTS_W Ptr<FastBilateralSolverFilter> createFastBilateralSolverFilter(InputArray guide, double sigma_spatial, double sigma_luma, double sigma_chroma, int num_iter = 25, double max_tol = 1e-5);
+
+
+/** @brief Simple one-line Fast Bilateral Solver filter call. If you have multiple images to filter with the same
+guide then use FastBilateralSolverFilter interface to avoid extra computations.
+
+@param guide image serving as guide for filtering. It should have 8-bit depth and either 1 or 3 channels.
+
+@param src source image for filtering with unsigned 8-bit or signed 16-bit or floating-point 32-bit depth and up to 4 channels.
+
+@param confidence confidence image with unsigned 8-bit or signed 16-bit or floating-point 32-bit confidence and 1 channel.
+
+@param dst destination image.
+
+@param sigma_spatial parameter, that is similar to spatial space sigma in bilateralFilter.
+
+@param sigma_luma parameter, that is similar to luma space sigma in bilateralFilter.
+
+@param sigma_chroma parameter, that is similar to chroma space sigma in bilateralFilter.
+
+@param num_iter number of iterations used for solving, 25 is usually enough.
+
+@param max_tol solving tolerance used for solving, 25 is usually enough.
+*/
+CV_EXPORTS_W void fastBilateralSolverFilter(InputArray guide, InputArray src, InputArray confidence, OutputArray dst, double sigma_spatial = 8, double sigma_luma = 8, double sigma_chroma = 8, int num_iter = 25, double max_tol = 1e-5);
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
 
 /** @brief Interface for implementations of Fast Global Smoother filter.
 

modules/ximgproc/samples/disparity_filtering.cpp

@@ -3,7 +3,7 @@
 #include "opencv2/imgcodecs.hpp"
 #include "opencv2/highgui.hpp"
 #include "opencv2/core/utility.hpp"
-#include "opencv2/ximgproc/disparity_filter.hpp"
+#include "opencv2/ximgproc.hpp"
 #include <iostream>
 #include <string>
 
@@ -21,7 +21,7 @@ const String keys =
     "{dst_path       |None              | optional path to save the resulting filtered disparity map        }"
     "{dst_raw_path   |None              | optional path to save raw disparity map before filtering          }"
     "{algorithm      |bm                | stereo matching method (bm or sgbm)                               }"
-    "{filter         |wls_conf          | used post-filtering (wls_conf or wls_no_conf)                     }"
+    "{filter         |wls_conf          | used post-filtering (wls_conf or wls_no_conf or fbs_conf)                     }"
     "{no-display     |                  | don't display results                                             }"
     "{no-downscale   |                  | force stereo matching on full-sized views to improve quality      }"
     "{dst_conf_path  |None              | optional path to save the confidence map used in filtering        }"
@@ -107,14 +107,15 @@ int main(int argc, char** argv)
         }
     }
 
-    Mat left_for_matcher, right_for_matcher;
+    Mat left_for_matcher, right_for_matcher, guide;
     Mat left_disp,right_disp;
-    Mat filtered_disp;
+    Mat filtered_disp,solved_disp;
     Mat conf_map = Mat(left.rows,left.cols,CV_8U);
     conf_map = Scalar(255);
     Rect ROI;
     Ptr<DisparityWLSFilter> wls_filter;
     double matching_time, filtering_time;
+    double solving_time = 0;
     if(max_disp<=0 || max_disp%16!=0)
     {
         cout<<"Incorrect max_disparity value: it should be positive and divisible by 16";
@@ -125,6 +126,7 @@ int main(int argc, char** argv)
         cout<<"Incorrect window_size value: it should be positive and odd";
         return -1;
     }
+
     if(filter=="wls_conf") // filtering with confidence (significantly better quality than wls_no_conf)
     {
         if(!no_downscale)
@@ -201,6 +203,99 @@ int main(int argc, char** argv)
             ROI = Rect(ROI.x*2,ROI.y*2,ROI.width*2,ROI.height*2);
         }
     }
+    else if(filter=="fbs_conf") // filtering with confidence (significantly better quality than wls_no_conf)
+    {
+        if(!no_downscale)
+        {
+            // downscale the views to speed-up the matching stage, as we will need to compute both left
+            // and right disparity maps for confidence map computation
+            //! [downscale_wls]
+            max_disp/=2;
+            if(max_disp%16!=0)
+                max_disp += 16-(max_disp%16);
+            resize(left ,left_for_matcher ,Size(),0.5,0.5);
+            resize(right,right_for_matcher,Size(),0.5,0.5);
+            //! [downscale_wls]
+        }
+        else
+        {
+            left_for_matcher  = left.clone();
+            right_for_matcher = right.clone();
+        }
+        guide = left_for_matcher.clone();
+
+        if(algo=="bm")
+        {
+            //! [matching_wls]
+            Ptr<StereoBM> left_matcher = StereoBM::create(max_disp,wsize);
+            wls_filter = createDisparityWLSFilter(left_matcher);
+            Ptr<StereoMatcher> right_matcher = createRightMatcher(left_matcher);
+
+            cvtColor(left_for_matcher,  left_for_matcher,  COLOR_BGR2GRAY);
+            cvtColor(right_for_matcher, right_for_matcher, COLOR_BGR2GRAY);
+
+            matching_time = (double)getTickCount();
+            left_matcher-> compute(left_for_matcher, right_for_matcher,left_disp);
+            right_matcher->compute(right_for_matcher,left_for_matcher, right_disp);
+            matching_time = ((double)getTickCount() - matching_time)/getTickFrequency();
+            //! [matching_wls]
+        }
+        else if(algo=="sgbm")
+        {
+            Ptr<StereoSGBM> left_matcher  = StereoSGBM::create(0,max_disp,wsize);
+            left_matcher->setP1(24*wsize*wsize);
+            left_matcher->setP2(96*wsize*wsize);
+            left_matcher->setPreFilterCap(63);
+            left_matcher->setMode(StereoSGBM::MODE_SGBM_3WAY);
+            wls_filter = createDisparityWLSFilter(left_matcher);
+            Ptr<StereoMatcher> right_matcher = createRightMatcher(left_matcher);
+
+            matching_time = (double)getTickCount();
+            left_matcher-> compute(left_for_matcher, right_for_matcher,left_disp);
+            right_matcher->compute(right_for_matcher,left_for_matcher, right_disp);
+            matching_time = ((double)getTickCount() - matching_time)/getTickFrequency();
+        }
+        else
+        {
+            cout<<"Unsupported algorithm";
+            return -1;
+        }
+
+        //! [filtering_wls]
+        wls_filter->setLambda(lambda);
+        wls_filter->setSigmaColor(sigma);
+        filtering_time = (double)getTickCount();
+        wls_filter->filter(left_disp,left,filtered_disp,right_disp);
+        filtering_time = ((double)getTickCount() - filtering_time)/getTickFrequency();
+        //! [filtering_wls]
+
+        conf_map = wls_filter->getConfidenceMap();
+
+        Mat left_disp_resized;
+        resize(left_disp,left_disp_resized,left.size());
+
+        // Get the ROI that was used in the last filter call:
+        ROI = wls_filter->getROI();
+        if(!no_downscale)
+        {
+            // upscale raw disparity and ROI back for a proper comparison:
+            resize(left_disp,left_disp,Size(),2.0,2.0);
+            left_disp = left_disp*2.0;
+            ROI = Rect(ROI.x*2,ROI.y*2,ROI.width*2,ROI.height*2);
+        }
+
+#ifdef HAVE_EIGEN
+        //! [filtering_fbs]
+        solving_time = (double)getTickCount();
+        // wls_filter->filter(left_disp,left,filtered_disp,right_disp);
+        // fastBilateralSolverFilter(left, filtered_disp, conf_map, solved_disp, 16.0, 16.0, 16.0);
+        fastBilateralSolverFilter(left, left_disp_resized, conf_map, solved_disp, 16.0, 16.0, 16.0);
+        solving_time = ((double)getTickCount() - solving_time)/getTickFrequency();
+	      solved_disp.convertTo(solved_disp, CV_8UC1);
+        cv::equalizeHist(solved_disp, solved_disp);
+        //! [filtering_fbs]
+#endif
+    }
     else if(filter=="wls_no_conf")
     {
         /* There is no convenience function for the case of filtering with no confidence, so we
@@ -263,6 +358,7 @@ int main(int argc, char** argv)
     cout.precision(2);
     cout<<"Matching time:  "<<matching_time<<"s"<<endl;
     cout<<"Filtering time: "<<filtering_time<<"s"<<endl;
+    cout<<"solving time: "<<solving_time<<"s"<<endl;
     cout<<endl;
 
     double MSE_before,percent_bad_before,MSE_after,percent_bad_after;
@@ -323,6 +419,33 @@ int main(int argc, char** argv)
         getDisparityVis(filtered_disp,filtered_disp_vis,vis_mult);
         namedWindow("filtered disparity", WINDOW_AUTOSIZE);
         imshow("filtered disparity", filtered_disp_vis);
+
+        if(!solved_disp.empty())
+        {
+            namedWindow("solved disparity", WINDOW_AUTOSIZE);
+            imshow("solved disparity", solved_disp);
+
+#define ENABLE_DOMAIN_TRANSFORM_FILTER
+#ifdef ENABLE_DOMAIN_TRANSFORM_FILTER
+	const float property_dt_sigmaSpatial = 40.0f;
+	const float property_dt_sigmaColor = 220.0f;
+	const int property_dt_numIters = 3;
+	cv::Mat final_disparty_dtfiltered_image;
+	cv::ximgproc::dtFilter(left,
+		solved_disp, final_disparty_dtfiltered_image,
+		property_dt_sigmaSpatial, property_dt_sigmaColor,
+		cv::ximgproc::DTF_RF,
+		property_dt_numIters);
+
+	// display disparity image
+	cv::Mat adjmap_dt;
+	final_disparty_dtfiltered_image.convertTo(adjmap_dt, CV_8UC1);
+		// 255.0f / 255.0f, 0.0f);
+	cv::imshow("disparity image + domain transform", adjmap_dt);
+#endif
+
+        }
+
         waitKey();
         //! [visualization]
     }

modules/ximgproc/src/fgs_filter.cpp

@@ -486,7 +486,6 @@ void FastGlobalSmootherFilterImpl::VerticalPass_ParBody::operator()(const Range&
     int start = std::min(range.start * stripe_sz, w);
     int end   = std::min(range.end   * stripe_sz, w);
 
-    //float lambda = fgs->lambda;
     WorkType denom;
     WorkType *Cvert_row, *Cvert_row_prev;
     WorkType *interD_row, *interD_row_prev, *cur_row, *cur_row_prev, *cur_row_next;
@@ -677,13 +676,11 @@ void FastGlobalSmootherFilterImpl::ComputeLUT_ParBody::operator()(const Range& r
 
 ////////////////////////////////////////////////////////////////////////////////////////////////
 
-CV_EXPORTS_W
 Ptr<FastGlobalSmootherFilter> createFastGlobalSmootherFilter(InputArray guide, double lambda, double sigma_color, double lambda_attenuation, int num_iter)
 {
     return Ptr<FastGlobalSmootherFilter>(FastGlobalSmootherFilterImpl::create(guide, lambda, sigma_color, num_iter, lambda_attenuation));
 }
 
-CV_EXPORTS_W
 void fastGlobalSmootherFilter(InputArray guide, InputArray src, OutputArray dst, double lambda, double sigma_color, double lambda_attenuation, int num_iter)
 {
     Ptr<FastGlobalSmootherFilter> fgs = createFastGlobalSmootherFilter(guide, lambda, sigma_color, lambda_attenuation, num_iter);

modules/ximgproc/test/test_fbs_filter.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"
+
+#ifdef HAVE_EIGEN
+
+namespace opencv_test { namespace {
+
+using namespace std;
+using namespace cv;
+using namespace cv::ximgproc;
+
+static string getDataDir()
+{
+    return cvtest::TS::ptr()->get_data_path();
+}
+
+CV_ENUM(SrcTypes, CV_8UC1, CV_8UC3, CV_8UC4, CV_16SC1, CV_16SC3, CV_32FC1);
+CV_ENUM(GuideTypes, CV_8UC1, CV_8UC3)
+typedef tuple<Size, SrcTypes, GuideTypes> FBSParams;
+typedef TestWithParam<FBSParams> FastBilateralSolverTest;
+
+TEST(FastBilateralSolverTest, SplatSurfaceAccuracy)
+{
+    RNG rnd(0);
+    int chanLut[] = {1,3,4};
+
+    for (int i = 0; i < 5; i++)
+    {
+        Size sz(rnd.uniform(512, 1024), rnd.uniform(512, 1024));
+
+        int guideCn = rnd.uniform(0, 2); // 1 or 3 channels
+        Mat guide(sz, CV_MAKE_TYPE(CV_8U, chanLut[guideCn]));
+        randu(guide, 0, 255);
+
+        Scalar surfaceValue;
+        int srcCn = rnd.uniform(0, 3); // 1, 3 or 4 channels
+        rnd.fill(surfaceValue, RNG::UNIFORM, 0, 255);
+        Mat src(sz, CV_MAKE_TYPE(CV_16S, chanLut[srcCn]), surfaceValue);
+        Mat confidence(sz, CV_MAKE_TYPE(CV_8U, 1), 255);
+
+        double sigma_spatial = rnd.uniform(4.0, 40.0);
+        double sigma_luma = rnd.uniform(4.0, 40.0);
+        double sigma_chroma  = rnd.uniform(4.0, 40.0);
+
+        Mat res;
+        fastBilateralSolverFilter(guide, src, confidence, res, sigma_spatial, sigma_luma, sigma_chroma);
+
+        // 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);
+    }
+}
+
+TEST(FastBilateralSolverTest, ReferenceAccuracy)
+{
+    string dir = getDataDir() + "cv/edgefilter";
+
+    Mat src = imread(dir + "/kodim23.png");
+    Mat ref = imread(dir + "/fgs/kodim23_lambda=1000_sigma=10.png");
+    Mat confidence(src.size(), CV_MAKE_TYPE(CV_8U, 1), 255);
+
+    ASSERT_FALSE(src.empty());
+    ASSERT_FALSE(ref.empty());
+
+    Mat res;
+    fastBilateralSolverFilter(src,src,confidence,res, 16.0, 16.0, 16.0);
+
+    double totalMaxError = 1.0/64.0*src.total()*src.channels();
+
+    EXPECT_LE(cvtest::norm(res, ref, NORM_L2), totalMaxError);
+    EXPECT_LE(cvtest::norm(res, ref, NORM_INF), 100);
+}
+
+INSTANTIATE_TEST_CASE_P(FullSet, FastBilateralSolverTest,Combine(Values(szODD, szQVGA), SrcTypes::all(), GuideTypes::all()));
+
+}
+}
+
+#endif //HAVE_EIGEN