Merge pull request #1877 from ytyytyyt:master

* fbs disparity filtering fix

* fbs fix vol II

* fbs doc update

* fix unused variables warnings in disparity_filtering.cpp

* trailing whitespaces removed

modules/ximgproc/doc/ximgproc.bib

@@ -311,9 +311,11 @@
   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},
+@inproceedings{BarronPoole2016,
+author = {Jonathan T Barron and Ben Poole},
+title={The Fast Bilateral Solver},
+booktitle={European Conference on Computer Vision (ECCV)},
+year={2016},
+publisher={Springer International Publishing},
+pages={617--632},
 }

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

@@ -384,9 +384,9 @@ 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 src source image for filtering with unsigned 8-bit or signed 16-bit or floating-point 32-bit depth and up to 3 channels.
 
-    @param confidence confidence image with unsigned 8-bit or signed 16-bit or floating-point 32-bit confidence and 1 channel.
+    @param confidence confidence image with unsigned 8-bit or floating-point 32-bit confidence and 1 channel.
 
     @param dst destination image.
     */
@@ -405,9 +405,10 @@ public:
 
 @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.
+@param max_tol solving tolerance used for solving.
 
 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);
 
@@ -419,7 +420,7 @@ guide then use FastBilateralSolverFilter interface to avoid extra computations.
 
 @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 confidence confidence image with unsigned 8-bit or floating-point 32-bit confidence and 1 channel.
 
 @param dst destination image.
 
@@ -431,7 +432,9 @@ guide then use FastBilateralSolverFilter interface to avoid extra computations.
 
 @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.
+@param max_tol solving tolerance used for solving.
+
+@note Confidence images with CV_8U depth are expected to in [0, 255] and CV_32F in [0, 1] range.
 */
 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);
 //////////////////////////////////////////////////////////////////////////

modules/ximgproc/samples/disparity_filtering.cpp

@@ -21,15 +21,18 @@ 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 or fbs_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        }"
     "{vis_mult       |1.0               | coefficient used to scale disparity map visualizations            }"
     "{max_disparity  |160               | parameter of stereo matching                                      }"
     "{window_size    |-1                | parameter of stereo matching                                      }"
-    "{wls_lambda     |8000.0            | parameter of post-filtering                                       }"
-    "{wls_sigma      |1.5               | parameter of post-filtering                                       }"
+    "{wls_lambda     |8000.0            | parameter of wls post-filtering                                   }"
+    "{wls_sigma      |1.5               | parameter of wls post-filtering                                   }"
+    "{fbs_spatial    |16.0              | parameter of fbs post-filtering                                   }"
+    "{fbs_luma       |8.0               | parameter of fbs post-filtering                                   }"
+    "{fbs_chroma     |8.0               | parameter of fbs post-filtering                                   }"
     ;
 
 int main(int argc, char** argv)
@@ -56,6 +59,9 @@ int main(int argc, char** argv)
     int max_disp = parser.get<int>("max_disparity");
     double lambda = parser.get<double>("wls_lambda");
     double sigma  = parser.get<double>("wls_sigma");
+    double fbs_spatial = parser.get<double>("fbs_spatial");
+    double fbs_luma = parser.get<double>("fbs_luma");
+    double fbs_chroma = parser.get<double>("fbs_chroma");
     double vis_mult = parser.get<double>("vis_mult");
 
     int wsize;
@@ -107,9 +113,9 @@ int main(int argc, char** argv)
         }
     }
 
-    Mat left_for_matcher, right_for_matcher, guide;
+    Mat left_for_matcher, right_for_matcher;
     Mat left_disp,right_disp;
-    Mat filtered_disp,solved_disp;
+    Mat filtered_disp,solved_disp,solved_filtered_disp;
     Mat conf_map = Mat(left.rows,left.cols,CV_8U);
     conf_map = Scalar(255);
     Rect ROI;
@@ -203,7 +209,7 @@ 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)
+    else if(filter=="fbs_conf") // filtering with fbs and confidence using also wls pre-processing
     {
         if(!no_downscale)
         {
@@ -222,7 +228,6 @@ int main(int argc, char** argv)
             left_for_matcher  = left.clone();
             right_for_matcher = right.clone();
         }
-        guide = left_for_matcher.clone();
 
         if(algo=="bm")
         {
@@ -281,20 +286,26 @@ int main(int argc, char** argv)
             // 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;
+            left_disp_resized = left_disp_resized*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);
+        fastBilateralSolverFilter(left, left_disp_resized, conf_map/255.0f, solved_disp, fbs_spatial, fbs_luma, fbs_chroma);
         solving_time = ((double)getTickCount() - solving_time)/getTickFrequency();
-	      solved_disp.convertTo(solved_disp, CV_8UC1);
-        cv::equalizeHist(solved_disp, solved_disp);
         //! [filtering_fbs]
+
+        //! [filtering_wls2fbs]
+        fastBilateralSolverFilter(left, filtered_disp, conf_map/255.0f, solved_filtered_disp, fbs_spatial, fbs_luma, fbs_chroma);
+        //! [filtering_wls2fbs]
+#else
+        (void)fbs_spatial;
+        (void)fbs_luma;
+        (void)fbs_chroma;
 #endif
+
     }
     else if(filter=="wls_no_conf")
     {
@@ -358,7 +369,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<<"Solving time: "<<solving_time<<"s"<<endl;
     cout<<endl;
 
     double MSE_before,percent_bad_before,MSE_after,percent_bad_after;
@@ -422,31 +433,23 @@ int main(int argc, char** argv)
 
         if(!solved_disp.empty())
         {
+            Mat solved_disp_vis;
+            getDisparityVis(solved_disp,solved_disp_vis,vis_mult);
             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
+            imshow("solved disparity", solved_disp_vis);
 
+            Mat solved_filtered_disp_vis;
+            getDisparityVis(solved_filtered_disp,solved_filtered_disp_vis,vis_mult);
+            namedWindow("solved wls disparity", WINDOW_AUTOSIZE);
+            imshow("solved wls disparity", solved_filtered_disp_vis);
         }
 
-        waitKey();
+        while(1)
+        {
+            char key = (char)waitKey();
+            if( key == 27 || key == 'q' || key == 'Q') // 'ESC'
+                break;
+        }
         //! [visualization]
     }
 

modules/ximgproc/src/fbs_filter.cpp

@@ -66,6 +66,8 @@
     typedef std::unordered_map<long long /* hash */, int /* vert id */>  mapId;
 #endif
 
+#define EPS 1e-43f
+
 namespace cv
 {
 namespace ximgproc
@@ -88,7 +90,7 @@ namespace ximgproc
         {
 
             CV_Assert(!src.empty() && (src.depth() == CV_8U || src.depth() == CV_16S || src.depth() == CV_32F) && src.channels()<=4);
-            CV_Assert(!confidence.empty() && (confidence.depth() == CV_8U || confidence.depth() == CV_16S || confidence.depth() == CV_32F) && confidence.channels()==1);
+            CV_Assert(!confidence.empty() && (confidence.depth() == CV_8U || confidence.depth() == CV_32F) && confidence.channels()==1);
             if (src.rows() != rows || src.cols() != cols)
             {
                 CV_Error(Error::StsBadSize, "Size of the filtered image must be equal to the size of the guide image");
@@ -108,14 +110,10 @@ namespace ximgproc
                 split(src,src_channels);
 
             Mat conf = confidence.getMat();
-            if(conf.depth() != CV_8UC1)
-                conf.convertTo(conf, CV_8UC1);
 
             for(int i=0;i<src.channels();i++)
             {
                 Mat cur_res = src_channels[i].clone();
-                if(src.depth() != CV_8UC1)
-                    cur_res.convertTo(cur_res, CV_8UC1);
 
                 solve(cur_res,conf,cur_res);
                 cur_res.convertTo(cur_res, src.type());
@@ -167,7 +165,6 @@ namespace ximgproc
         Eigen::SparseMatrix<float, Eigen::ColMajor> S;
         Eigen::SparseMatrix<float, Eigen::ColMajor> Dn;
         Eigen::SparseMatrix<float, Eigen::ColMajor> Dm;
-        cv::Mat guide;
 
         struct grid_params
         {
@@ -206,7 +203,6 @@ namespace ximgproc
 
     void FastBilateralSolverFilterImpl::init(cv::Mat& reference, double sigma_spatial, double sigma_luma, double sigma_chroma, int num_iter, double max_tol)
     {
-        guide = reference.clone();
 
         bs_param.cg_maxiter = num_iter;
         bs_param.cg_tol = max_tol;
@@ -479,45 +475,69 @@ namespace ximgproc
         Eigen::VectorXf y1(nvertices);
         Eigen::VectorXf w_splat(nvertices);
 
-        cv::Mat x;
-        cv::Mat w;
-        cv::Mat xw;
-        cv::Mat filtered_xw;
-        cv::Mat filtered_w;
-        cv::Mat filtered_disp;
-        float fgs_colorSigma = 1.5;
-        float fgs_spatialSigma = 2000;
-	      target.convertTo(x, CV_32FC1, 1.0f/255.0f);
-	      confidence.convertTo(w, CV_32FC1);
-        xw = x.mul(w);
-        cv::ximgproc::fastGlobalSmootherFilter(guide, xw, filtered_xw, fgs_spatialSigma, fgs_colorSigma);
-        cv::ximgproc::fastGlobalSmootherFilter(guide, w, filtered_w, fgs_spatialSigma, fgs_colorSigma);
-        cv::divide(filtered_xw, filtered_w, filtered_disp, 1.0f, CV_32FC1);
+        Eigen::VectorXf x(npixels);
+        Eigen::VectorXf w(npixels);
 
+        if(target.depth() == CV_16S)
+        {
+            const int16_t *pft = reinterpret_cast<const int16_t*>(target.data);
+            for (int i = 0; i < npixels; i++)
+            {
+                x(i) = (cv::saturate_cast<float>(pft[i])+32768.0f)/65535.0f;
+            }
+        }
+        else if(target.depth() == CV_8U)
+        {
+            const uchar *pft = reinterpret_cast<const uchar*>(target.data);
+            for (int i = 0; i < npixels; i++)
+            {
+                x(i) = cv::saturate_cast<float>(pft[i])/255.0f;
+            }
+        }
+        else if(confidence.depth() == CV_32F)
+        {
+            const float *pft = reinterpret_cast<const float*>(target.data);
+            for (int i = 0; i < npixels; i++)
+            {
+                x(i) = pft[i];
+            }
+        }
 
-        //construct A
-        w_splat.setZero();
-        const float *pfw = reinterpret_cast<const float*>(w.data);
-        for (int i = 0; i < int(splat_idx.size()); i++)
+        if(confidence.depth() == CV_8U)
         {
-            w_splat(splat_idx[i]) += pfw[i]/255.0f;
+            const uchar *pfc = reinterpret_cast<const uchar*>(confidence.data);
+            for (int i = 0; i < npixels; i++)
+            {
+                w(i) = cv::saturate_cast<float>(pfc[i])/255.0f;
+            }
         }
+        else if(confidence.depth() == CV_32F)
+        {
+            const float *pfc = reinterpret_cast<const float*>(confidence.data);
+            for (int i = 0; i < npixels; i++)
+            {
+                w(i) = pfc[i];
+            }
+        }
+
+        //construct A
+        Splat(w,w_splat);
+
         diagonal(w_splat,A_data);
         A = bs_param.lam * (Dm - Dn * (blurs*Dn)) + A_data ;
 
         //construct b
         b.setZero();
-        const float *pfx = reinterpret_cast<const float*>(filtered_disp.data);
         for (int i = 0; i < int(splat_idx.size()); i++)
         {
-            b(splat_idx[i]) += pfx[i]*pfw[i]/255.0f;
+            b(splat_idx[i]) += x(i) * w(i);
         }
 
         //construct guess for y
         y0.setZero();
         for (int i = 0; i < int(splat_idx.size()); i++)
         {
-            y0(splat_idx[i]) += pfx[i];
+            y0(splat_idx[i]) += x(i);
         }
         y1.setZero();
         for (int i = 0; i < int(splat_idx.size()); i++)
@@ -541,10 +561,29 @@ namespace ximgproc
         std::cout << "estimated error: " << cg.error()      << std::endl;
 
         //slice
-        uchar *pftar = (uchar*)(output.data);
-        for (int i = 0; i < int(splat_idx.size()); i++)
+        if(target.depth() == CV_16S)
+        {
+            int16_t *pftar = (int16_t*) output.data;
+            for (int i = 0; i < int(splat_idx.size()); i++)
+            {
+                pftar[i] = cv::saturate_cast<ushort>(y(splat_idx[i]) * 65535.0f - 32768.0f);
+            }
+        }
+        else if (target.depth() == CV_8U)
         {
-            pftar[i] = uchar(y(splat_idx[i]) * 255.0f);
+            uchar *pftar = (uchar*) output.data;
+            for (int i = 0; i < int(splat_idx.size()); i++)
+            {
+                pftar[i] = cv::saturate_cast<uchar>(y(splat_idx[i]) * 255.0f);
+            }
+        }
+        else
+        {
+            float *pftar = (float*)(output.data);
+            for (int i = 0; i < int(splat_idx.size()); i++)
+            {
+                pftar[i] = y(splat_idx[i]);
+            }
         }