added DisparityBilateralFilter to gpu module

modules/gpu/include/opencv2/gpu/gpu.hpp

@@ -391,6 +391,7 @@ namespace cv
         ////////////////////////// StereoBeliefPropagation ///////////////////////////
         // "Efficient Belief Propagation for Early Vision" 
         // P.Felzenszwalb
+
         class CV_EXPORTS StereoBeliefPropagation
         {
         public:
@@ -504,6 +505,45 @@ namespace cv
 
             GpuMat out;
         };
+
+        /////////////////////////// DisparityBilateralFilter ///////////////////////////
+        // Disparity map refinement using joint bilateral filtering given a single color image.
+        // Qingxiong Yang, Liang Wang†, Narendra Ahuja         
+        // http://vision.ai.uiuc.edu/~qyang6/
+
+        class CV_EXPORTS DisparityBilateralFilter
+        {
+        public:
+            enum { DEFAULT_NDISP  = 64 };
+            enum { DEFAULT_RADIUS = 3 };
+            enum { DEFAULT_ITERS  = 1 };
+
+            //! the default constructor
+            explicit DisparityBilateralFilter(int ndisp = DEFAULT_NDISP, int radius = DEFAULT_RADIUS, int iters = DEFAULT_ITERS);
+
+            //! the full constructor taking the number of disparities, filter radius,
+            //! number of iterations, truncation of data continuity, truncation of disparity continuity
+            //! and filter range sigma
+            DisparityBilateralFilter(int ndisp, int radius, int iters, float edge_threshold, float max_disc_threshold, float sigma_range);
+
+            //! the disparity map refinement operator. Refine disparity map using joint bilateral filtering given a single color image.
+            //! disparity must have CV_8U or CV_16S type, image must have CV_8UC1 or CV_8UC3 type.
+            void operator()(const GpuMat& disparity, const GpuMat& image, GpuMat& dst);
+
+            //! Acync version
+            void operator()(const GpuMat& disparity, const GpuMat& image, GpuMat& dst, Stream& stream);
+
+            int ndisp;
+            int radius;
+            int iters;
+
+            float edge_threshold;
+            float max_disc_threshold;
+            float sigma_range;
+        private:
+            std::vector<float> table_color;
+            GpuMat table_space;
+        };
     }
 
     //! Speckle filtering - filters small connected components on diparity image.

modules/gpu/src/bilateral_filter.cpp

+/*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, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * 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 GpuMaterials provided with the distribution.
+//
+//   * The name of the copyright holders 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"
+
+using namespace cv;
+using namespace cv::gpu;
+using namespace std;
+
+#if !defined (HAVE_CUDA)
+
+cv::gpu::DisparityBilateralFilter::DisparityBilateralFilter(int, int, int) { throw_nogpu(); }
+cv::gpu::DisparityBilateralFilter::DisparityBilateralFilter(int, int, int, float, float, float) { throw_nogpu(); }
+
+void cv::gpu::DisparityBilateralFilter::operator()(const GpuMat&, const GpuMat&, GpuMat&) { throw_nogpu(); }
+void cv::gpu::DisparityBilateralFilter::operator()(const GpuMat&, const GpuMat&, GpuMat&, Stream&) { throw_nogpu(); }
+
+#else /* !defined (HAVE_CUDA) */
+
+namespace cv { namespace gpu { namespace bf 
+{
+    void calc_space_weighted_filter_gpu(const DevMem2Df& table_space, int half, float dist_space, cudaStream_t stream);
+    void load_constants(float* table_color, const DevMem2Df& table_space, int ndisp, int radius, short edge_disc, short max_disc);
+
+    void bilateral_filter_gpu(const DevMem2D& disp, const DevMem2D& img, int channels, int iters, cudaStream_t stream);
+    void bilateral_filter_gpu(const DevMem2D_<short>& disp, const DevMem2D& img, int channels, int iters, cudaStream_t stream);
+}}}
+
+namespace
+{
+    const float DEFAULT_EDGE_THRESHOLD = 0.1f;
+    const float DEFAULT_MAX_DISC_THRESHOLD = 0.2f;
+    const float DEFAULT_SIGMA_RANGE = 10.0f;
+
+    inline void calc_color_weighted_table(vector<float>& table_color, float sigma_range, int len)
+    {
+        float* color_table_x;
+        
+        table_color.resize(len);
+        color_table_x = &table_color[0];
+
+        for(int y = 0; y < len; y++) 
+            table_color[y] = static_cast<float>(std::exp(-double(y * y) / (2 * sigma_range * sigma_range)));
+    }
+
+    inline void calc_space_weighted_filter(GpuMat& table_space, int win_size, float dist_space, cudaStream_t stream)
+    {
+        int half = (win_size >> 1);
+        table_space.create(half + 1, half + 1, CV_32F);
+
+        bf::calc_space_weighted_filter_gpu(table_space, half, dist_space, stream);
+    }
+
+    template <typename T>
+    void bilateral_filter_operator(DisparityBilateralFilter& rthis, vector<float>& table_color, GpuMat& table_space, 
+                                   const GpuMat& disp, const GpuMat& img, GpuMat& dst, cudaStream_t stream)
+    {
+        calc_color_weighted_table(table_color, rthis.sigma_range, 255);
+        calc_space_weighted_filter(table_space, rthis.radius * 2 + 1, rthis.radius + 1.0f, stream);
+
+        short edge_disc = max<short>(short(1), short(rthis.ndisp * rthis.edge_threshold + 0.5));
+        short max_disc = short(rthis.ndisp * rthis.max_disc_threshold + 0.5);
+
+        bf::load_constants(&table_color[0], table_space, rthis.ndisp, rthis.radius, edge_disc, max_disc);
+
+        if (&dst != &disp)
+            disp.copyTo(dst);
+
+        bf::bilateral_filter_gpu((DevMem2D_<T>)dst, img, img.channels(), rthis.iters, stream);
+    }
+
+    typedef void (*bilateral_filter_operator_t)(DisparityBilateralFilter& rthis, vector<float>& table_color, GpuMat& table_space, 
+                                   const GpuMat& disp, const GpuMat& img, GpuMat& dst, cudaStream_t stream);
+    
+    const bilateral_filter_operator_t operators[] = 
+        {bilateral_filter_operator<unsigned char>, 0, 0, bilateral_filter_operator<short>, 0, 0, 0, 0};
+}
+
+cv::gpu::DisparityBilateralFilter::DisparityBilateralFilter(int ndisp_, int radius_, int iters_)
+    : ndisp(ndisp_), radius(radius_), iters(iters_), edge_threshold(DEFAULT_EDGE_THRESHOLD), max_disc_threshold(DEFAULT_MAX_DISC_THRESHOLD),
+      sigma_range(DEFAULT_SIGMA_RANGE)
+{
+}
+
+cv::gpu::DisparityBilateralFilter::DisparityBilateralFilter(int ndisp_, int radius_, int iters_, float edge_threshold_, 
+                                                     float max_disc_threshold_, float sigma_range_)
+    : ndisp(ndisp_), radius(radius_), iters(iters_), edge_threshold(edge_threshold_), max_disc_threshold(max_disc_threshold_), 
+      sigma_range(sigma_range_)
+{
+}
+
+void cv::gpu::DisparityBilateralFilter::operator()(const GpuMat& disp, const GpuMat& img, GpuMat& dst)
+{
+    CV_DbgAssert(0 < ndisp && 0 < radius && 0 < iters);
+    CV_Assert(disp.rows == img.rows && disp.cols == img.cols && (disp.type() == CV_8U || disp.type() == CV_16S) && (img.type() == CV_8UC1 || img.type() == CV_8UC3));
+    operators[disp.type()](*this, table_color, table_space, disp, img, dst, 0);
+}
+
+void cv::gpu::DisparityBilateralFilter::operator()(const GpuMat& disp, const GpuMat& img, GpuMat& dst, Stream& stream)
+{
+    CV_DbgAssert(0 < ndisp && 0 < radius && 0 < iters);
+    CV_Assert(disp.rows == img.rows && disp.cols == img.cols && (disp.type() == CV_8U || disp.type() == CV_16S) && (img.type() == CV_8UC1 || img.type() == CV_8UC3));
+    operators[disp.type()](*this, table_color, table_space, disp, img, dst, StreamAccessor::getStream(stream));
+}
+
+#endif /* !defined (HAVE_CUDA) */

modules/gpu/src/cuda/bilateral_filter.cu

+/*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, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * 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 the copyright holders 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 "opencv2/gpu/devmem2d.hpp"
+#include "saturate_cast.hpp"
+#include "safe_call.hpp"
+
+using namespace cv::gpu;
+using namespace cv::gpu::impl;
+
+#ifndef FLT_MAX
+#define FLT_MAX 3.402823466e+30F
+#endif
+
+namespace bf_krnls
+{
+    __global__ void calc_space_weighted_filter(float* table_space, size_t step, int half, float dist_space)
+    {
+        int x = blockIdx.x * blockDim.x + threadIdx.x;
+        int y = blockIdx.y * blockDim.y + threadIdx.y;
+
+        if (y <= half && x <= half)
+            *(table_space + y * step + x) = expf(-sqrtf(float(y * y) + float(x * x)) / dist_space);
+    }
+}
+
+namespace cv { namespace gpu { namespace bf 
+{
+    void calc_space_weighted_filter_gpu(const DevMem2Df& table_space, int half, float dist_space, cudaStream_t stream)
+    {
+        dim3 threads(32, 8, 1);
+        dim3 grid(1, 1, 1);
+        grid.x = divUp(half + 1, threads.x);
+        grid.x = divUp(half + 1, threads.y);
+
+        bf_krnls::calc_space_weighted_filter<<<grid, threads, 0, stream>>>(table_space.ptr, table_space.step/sizeof(float), half, dist_space);
+
+        if (stream != 0)
+            cudaSafeCall( cudaThreadSynchronize() );
+    }
+}}}
+
+namespace bf_krnls
+{
+    __constant__ float* ctable_color;
+    __constant__ float* ctable_space;
+    __constant__ size_t ctable_space_step;
+
+    __constant__ int cndisp;
+    __constant__ int cradius;
+
+    __constant__ short cedge_disc;
+    __constant__ short cmax_disc;
+}
+
+namespace cv { namespace gpu { namespace bf 
+{
+    void load_constants(float* table_color, const DevMem2Df& table_space, int ndisp, int radius, short edge_disc, short max_disc)
+    {
+        cudaSafeCall( cudaMemcpyToSymbol(bf_krnls::ctable_color, &table_color, sizeof(table_color)) );
+        cudaSafeCall( cudaMemcpyToSymbol(bf_krnls::ctable_space, &table_space.ptr, sizeof(table_space.ptr)) );
+        size_t table_space_step = table_space.step / sizeof(float);
+        cudaSafeCall( cudaMemcpyToSymbol(bf_krnls::ctable_space_step, &table_space_step, sizeof(size_t)) );
+        
+        cudaSafeCall( cudaMemcpyToSymbol(bf_krnls::cndisp, &ndisp, sizeof(int)) );
+        cudaSafeCall( cudaMemcpyToSymbol(bf_krnls::cradius, &radius, sizeof(int)) );
+        
+        cudaSafeCall( cudaMemcpyToSymbol(bf_krnls::cedge_disc, &edge_disc, sizeof(short)) );
+        cudaSafeCall( cudaMemcpyToSymbol(bf_krnls::cmax_disc, &max_disc, sizeof(short)) );
+    }
+}}}
+
+namespace bf_krnls
+{
+    template <int channels>
+    struct DistRgbMax
+    {
+        static __device__ uchar calc(const uchar* a, const uchar* b)
+        {
+            uchar x = abs(a[0] - b[0]);
+            uchar y = abs(a[1] - b[1]);
+            uchar z = abs(a[2] - b[2]);
+            return (max(max(x, y), z));
+        }
+    };
+
+    template <>
+    struct DistRgbMax<1>
+    {
+        static __device__ uchar calc(const uchar* a, const uchar* b)
+        {
+            return abs(a[0] - b[0]);
+        }
+    };
+
+    template <int channels, typename T>
+    __global__ void bilateral_filter(int t, T* disp, size_t disp_step, const uchar* img, size_t img_step, int h, int w)
+    {
+        const int y = blockIdx.y * blockDim.y + threadIdx.y;
+        const int x = ((blockIdx.x * blockDim.x + threadIdx.x) << 1) + ((y + t) & 1);
+
+        T dp[5];
+
+        if (y > 0 && y < h - 1 && x > 0 && x < w - 1)
+        {
+            dp[0] = *(disp + (y  ) * disp_step + x + 0);
+            dp[1] = *(disp + (y-1) * disp_step + x + 0);
+            dp[2] = *(disp + (y  ) * disp_step + x - 1);
+            dp[3] = *(disp + (y+1) * disp_step + x + 0);
+            dp[4] = *(disp + (y  ) * disp_step + x + 1);
+
+            if(abs(dp[1] - dp[0]) >= cedge_disc || abs(dp[2] - dp[0]) >= cedge_disc || abs(dp[3] - dp[0]) >= cedge_disc || abs(dp[4] - dp[0]) >= cedge_disc)            
+            {
+                const int ymin = max(0, y - cradius);
+                const int xmin = max(0, x - cradius);
+                const int ymax = min(h - 1, y + cradius);
+                const int xmax = min(w - 1, x + cradius);
+
+                float cost[] = {0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+                const uchar* ic = img + y * img_step + channels * x;
+
+                for(int yi = ymin; yi <= ymax; yi++)
+                {
+                    const T* disp_y = disp + yi * disp_step;
+
+                    for(int xi = xmin; xi <= xmax; xi++)
+                    {
+                        const uchar* in = img + yi * img_step + channels * xi;
+
+                        uchar dist_rgb = DistRgbMax<channels>::calc(in, ic);
+
+                        const float weight = ctable_color[dist_rgb] * (ctable_space + abs(y-yi)* ctable_space_step)[abs(x-xi)];
+
+                        const T disp_reg = disp_y[xi];
+
+                        cost[0] += min(cmax_disc, abs(disp_reg - dp[0])) * weight;
+                        cost[1] += min(cmax_disc, abs(disp_reg - dp[1])) * weight;
+                        cost[2] += min(cmax_disc, abs(disp_reg - dp[2])) * weight;
+                        cost[3] += min(cmax_disc, abs(disp_reg - dp[3])) * weight;
+                        cost[4] += min(cmax_disc, abs(disp_reg - dp[4])) * weight;
+                    }
+                }
+
+                float minimum = FLT_MAX;
+                int id = 0;
+
+                if (cost[0] < minimum)
+                {
+                    minimum = cost[0];
+                    id = 0;
+                }
+                if (cost[1] < minimum)
+                {
+                    minimum = cost[1];
+                    id = 1;
+                }
+                if (cost[2] < minimum)
+                {
+                    minimum = cost[2];
+                    id = 2;
+                }
+                if (cost[3] < minimum)
+                {
+                    minimum = cost[3];
+                    id = 3;
+                }
+                if (cost[4] < minimum)
+                {
+                    minimum = cost[4];
+                    id = 4;
+                }
+
+                *(disp + y * disp_step + x) = dp[id];
+            }
+        }
+    }
+}
+
+namespace cv { namespace gpu { namespace bf 
+{
+    template <typename T>     
+    void bilateral_filter_caller(const DevMem2D_<T>& disp, const DevMem2D& img, int channels, int iters, cudaStream_t stream)
+    {
+        dim3 threads(32, 8, 1);
+        dim3 grid(1, 1, 1);
+        grid.x = divUp(disp.cols, threads.x << 1);
+        grid.y = divUp(disp.rows, threads.y);
+
+        switch (channels)
+        {
+        case 1:
+            for (int i = 0; i < iters; ++i)
+            {
+                bf_krnls::bilateral_filter<1><<<grid, threads, 0, stream>>>(0, disp.ptr, disp.step/sizeof(T), img.ptr, img.step, disp.rows, disp.cols);
+                bf_krnls::bilateral_filter<1><<<grid, threads, 0, stream>>>(1, disp.ptr, disp.step/sizeof(T), img.ptr, img.step, disp.rows, disp.cols);
+            }
+            break;
+        case 3:
+            for (int i = 0; i < iters; ++i)
+            {
+                bf_krnls::bilateral_filter<3><<<grid, threads, 0, stream>>>(0, disp.ptr, disp.step/sizeof(T), img.ptr, img.step, disp.rows, disp.cols);
+                bf_krnls::bilateral_filter<3><<<grid, threads, 0, stream>>>(1, disp.ptr, disp.step/sizeof(T), img.ptr, img.step, disp.rows, disp.cols);
+            }
+            break;
+        default:
+            cv::gpu::error("Unsupported channels count", __FILE__, __LINE__);
+        }        
+
+        if (stream != 0)
+            cudaSafeCall( cudaThreadSynchronize() );
+    }
+
+    void bilateral_filter_gpu(const DevMem2D& disp, const DevMem2D& img, int channels, int iters, cudaStream_t stream)
+    {
+        bilateral_filter_caller(disp, img, channels, iters, stream);
+    }
+
+    void bilateral_filter_gpu(const DevMem2D_<short>& disp, const DevMem2D& img, int channels, int iters, cudaStream_t stream)
+    {
+        bilateral_filter_caller(disp, img, channels, iters, stream);
+    }
+}}}