From 1f04ea477f07192a20dd4c4d64857c65a2790a2b Mon Sep 17 00:00:00 2001 From: Vladislav Vinogradov Date: Thu, 19 Aug 2010 08:44:06 +0000 Subject: [PATCH] added DisparityBilateralFilter to gpu module --- modules/gpu/include/opencv2/gpu/gpu.hpp | 40 ++++ modules/gpu/src/bilateral_filter.cpp | 145 +++++++++++++ modules/gpu/src/cuda/bilateral_filter.cu | 262 +++++++++++++++++++++++ 3 files changed, 447 insertions(+) create mode 100644 modules/gpu/src/bilateral_filter.cpp create mode 100644 modules/gpu/src/cuda/bilateral_filter.cu diff --git a/modules/gpu/include/opencv2/gpu/gpu.hpp b/modules/gpu/include/opencv2/gpu/gpu.hpp index bf9abbaa0..e5c475f42 100644 --- a/modules/gpu/include/opencv2/gpu/gpu.hpp +++ b/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 table_color; + GpuMat table_space; + }; } //! Speckle filtering - filters small connected components on diparity image. diff --git a/modules/gpu/src/bilateral_filter.cpp b/modules/gpu/src/bilateral_filter.cpp new file mode 100644 index 000000000..ffe6a1e79 --- /dev/null +++ b/modules/gpu/src/bilateral_filter.cpp @@ -0,0 +1,145 @@ +/*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_& 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& 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(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 + void bilateral_filter_operator(DisparityBilateralFilter& rthis, vector& 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(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_)dst, img, img.channels(), rthis.iters, stream); + } + + typedef void (*bilateral_filter_operator_t)(DisparityBilateralFilter& rthis, vector& 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, 0, 0, bilateral_filter_operator, 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) */ diff --git a/modules/gpu/src/cuda/bilateral_filter.cu b/modules/gpu/src/cuda/bilateral_filter.cu new file mode 100644 index 000000000..f53d4c2b9 --- /dev/null +++ b/modules/gpu/src/cuda/bilateral_filter.cu @@ -0,0 +1,262 @@ +/*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<<>>(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 + 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 + __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::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 + void bilateral_filter_caller(const DevMem2D_& 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><<>>(0, disp.ptr, disp.step/sizeof(T), img.ptr, img.step, disp.rows, disp.cols); + bf_krnls::bilateral_filter<1><<>>(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><<>>(0, disp.ptr, disp.step/sizeof(T), img.ptr, img.step, disp.rows, disp.cols); + bf_krnls::bilateral_filter<3><<>>(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_& disp, const DevMem2D& img, int channels, int iters, cudaStream_t stream) + { + bilateral_filter_caller(disp, img, channels, iters, stream); + } +}}}