双边滤波CUDA优化——BilateralFilter CUDA

=======双边滤波概述=======

双边滤波（Bilateral filter）是一种可以保边去噪的滤波器。之所以可以达到此去噪效果，是因为滤波器是由两个函数构成。一个函数是由几何空间距离决定滤波器系数。另一个由像素差值决定滤波器系数。可以与其相比较的两个filter：高斯低通滤波器(http://en.wikipedia.org/wiki/Gaussian_filter)和α-截尾均值滤波器（去掉百分率为α的最小值和最大之后剩下像素的均值作为滤波器）。

=======双边滤波公式=======

=======双边滤波代码（CPU）=======

OpenCV源码：

/****************************************************************************************                                   Bilateral Filtering
\****************************************************************************************/

namespace cv
{

static void
bilateralFilter_8u( const Mat& src, Mat& dst, int d,
                    double sigma_color, double sigma_space,
                    int borderType )
{
    int cn = src.channels();
    int i, j, k, maxk, radius;
    Size size = src.size();

    CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) &&
        src.type() == dst.type() && src.size() == dst.size() &&
        src.data != dst.data );

    if( sigma_color <= 0 )
        sigma_color = 1;
    if( sigma_space <= 0 )
        sigma_space = 1;

    double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
    double gauss_space_coeff = -0.5/(sigma_space*sigma_space);

    if( d <= 0 )
        radius = cvRound(sigma_space*1.5);
    else
        radius = d/2;
    radius = MAX(radius, 1);
    d = radius*2 + 1;

    Mat temp;
    copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );

    vector<float> _color_weight(cn*256);
    vector<float> _space_weight(d*d);
    vector<int> _space_ofs(d*d);
    float* color_weight = &_color_weight[0];
    float* space_weight = &_space_weight[0];
    int* space_ofs = &_space_ofs[0];

    // initialize color-related bilateral filter coefficients
    for( i = 0; i < 256*cn; i++ )
        color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);

    // initialize space-related bilateral filter coefficients
    for( i = -radius, maxk = 0; i <= radius; i++ )
        for( j = -radius; j <= radius; j++ )
        {
            double r = std::sqrt((double)i*i + (double)j*j);
            if( r > radius )
                continue;
            space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
            space_ofs[maxk++] = (int)(i*temp.step + j*cn);
        }

    for( i = 0; i < size.height; i++ )
    {
        const uchar* sptr = temp.data + (i+radius)*temp.step + radius*cn;
        uchar* dptr = dst.data + i*dst.step;

        if( cn == 1 )
        {
            for( j = 0; j < size.width; j++ )
            {
                float sum = 0, wsum = 0;
                int val0 = sptr[j];
                for( k = 0; k < maxk; k++ )
                {
                    int val = sptr[j + space_ofs[k]];
                    float w = space_weight[k]*color_weight[std::abs(val - val0)];
                    sum += val*w;
                    wsum += w;
                }
                // overflow is not possible here => there is no need to use CV_CAST_8U
                dptr[j] = (uchar)cvRound(sum/wsum);
            }
        }
        else
        {
            assert( cn == 3 );
            for( j = 0; j < size.width*3; j += 3 )
            {
                float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
                int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
                for( k = 0; k < maxk; k++ )
                {
                    const uchar* sptr_k = sptr + j + space_ofs[k];
                    int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
                    float w = space_weight[k]*color_weight[std::abs(b - b0) +
                        std::abs(g - g0) + std::abs(r - r0)];
                    sum_b += b*w; sum_g += g*w; sum_r += r*w;
                    wsum += w;
                }
                wsum = 1.f/wsum;
                b0 = cvRound(sum_b*wsum);
                g0 = cvRound(sum_g*wsum);
                r0 = cvRound(sum_r*wsum);
                dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0;
            }
        }
    }
}

OpenCV 双边滤波调用

bilateralFilter(InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace,
                      int borderType=BORDER_DEFAULT )；

d 表示滤波时像素邻域直径，d为负时由 sigaColor计算得到；d>5时不能实时处理。

sigmaColor、sigmaSpace非别表示颜色空间和坐标空间的滤波系数sigma。可以简单的赋值为相同的值。<10时几乎没有效果；>150时为油画的效果。

borderType可以不指定。

OpenCV 双边滤波实验

用sigma为10，150,240,480 时效果如下：

      

=======双边滤波优化（CUDA）=======

在进行图像处理时，由于计算量大，常常无法到达实时的效果，因此需利用GPU处理，使用CUDA进行优化。尤其是图像滤波这种，(1) 并行度高，线程间耦合度低，每个像素的处理并不相互影响；(2) 像素传输量小，计算量大；特别适合CUDA进行计算。

CUDA BilateralFilter流程（可扩展至CUDA图像处理领域）

  1. 复制数据 Copy Data to Device

    • 在Device上开辟2维数据空间作为输入数据：

      • UINT *dImage = NULL; //original image
      • size_t pitch;
      • checkCudaErrors( cudaMallocPitch(&dImage, &pitch, sizeof(UINT)*width, height) );
    • 复制数据到显卡 Copy Data from Host to Device：
      • checkCudaErrors( cudaMemcpy2D(dImage, pitch, hImage, sizeof(UINT)*width, sizeof(UINT)*width,
        height, cudaMemcpyHostToDevice));
    • 在Device上开辟2维数据空间保存输出数据：
      • unsigned int *dResult;
      • size_t pitch;
      • checkCudaErrors( cudaMallocPitch((void **)&dResult, &pitch, width*sizeof(UINT), height)
        );
  2. 使用纹理储存器
    • 声明CUDA数组前，以结构体channelDesc描述CUDA数组中组件的数量和数据类型

      • cudaChannelFormatDesc desc = cudaCreateChannelDesc<uchar4>();
    • 声明纹理参照系:texture<Type,Dim,ReadMode> texRef;
      • texture<uchar4, 2, cudaReadModeElementType> rgbaTex; //以全局变量形式出现
    • 将纹理参照系绑定到一个CUDA数组
      • checkCudaErrors( cudaBindTexture2D(0, rgbaTex, dImage, desc, width, height, pitch));
    • 纹理拾取 // 书P65
      • tex1Dfetch()
      • tex1D(); tex2D(); tex3D();
  3. 使用常量储存器
    • 位于显存，但拥有缓存加速。空间较小（64K），保存频繁访问的只读数据，
    • 两种使用方法
      • 直接赋值：__constant__ float c_cuda[4] = {0,1,2,3};  __constant__ float c_num = 1;
      • 使用函数：__constant__ float color_weight[4*256];  checkCudaErrors(cudaMemcpyToSymbol(color_weight, color_gaussian, sizeof(float)*(4*256)));
    • 本程序
      • __constant__ float color_weight[4*256];
      • __constant__ float space_weight[1024];
      • 函数声明局部变量数组 float color_gaussian[4*256];  float space_gaussian[1024]; 使用定义域核和值域核赋值
      • 函数赋值：
        • checkCudaErrors(cudaMemcpyToSymbol(color_weight, color_gaussian, sizeof(float)*(4*256)));
        • checkCudaErrors(cudaMemcpyToSymbol(space_weight, space_gaussian, sizeof(float)*(1024)));
  4. CUDA处理Kernel
    • 并行设计：

      • dim3 blockSize(16, 16);
      • dim3 gridSize((width + 16 - 1) / 16, (height + 16 - 1) / 16);
      • d_bilateral_filter<<< gridSize, blockSize>>>(dDest, width, height, radius);
    • 纹理拾取
      • center = tex2D(rgbaTex, x, y);
    • 计算定义域和值域核的值
  5. 复制数据Copy Data to Host
    • checkCudaErrors( cudaMemcpy(hImage,dResult,(width * height)*sizeof(UINT),cudaMemcpyDeviceToHost));

CUDA源码：

BilateralKernel.cu

#include <helper_math.h>
#include <helper_functions.h>
#include <helper_cuda.h>       // CUDA device initialization helper functions
#include "cuda.h"
#include "cuda_runtime_api.h"
#include "device_launch_parameters.h"
#include "device_functions.h"

__constant__ float color_weight[4*256];
__constant__ float space_weight[1024];

UINT *dImage  = NULL;   //original image
UINT *dTemp   = NULL;   //temp array for iterations
size_t pitch;

texture<uchar4, 2, cudaReadModeElementType> rgbaTex;//声明纹理参照系

__device__ float colorLenGaussian(uchar4 a, uchar4 b)
{
	//若想达到漫画效果，就注释掉sqrt,使颜色距离变大
    uint mod = (uint)sqrt(((float)b.x - (float)a.x) * ((float)b.x - (float)a.x) +
				((float)b.y - (float)a.y) * ((float)b.y - (float)a.y) +
                ((float)b.z - (float)a.z) * ((float)b.z - (float)a.z) +
				((float)b.w - (float)a.w) * ((float)b.w - (float)a.w));

    return color_weight[mod];
}
__device__ uint rgbaFloatToInt(float4 rgba)
{
    rgba.x = __saturatef(fabs(rgba.x));   // clamp to [0.0, 1.0]
    rgba.y = __saturatef(fabs(rgba.y));
    rgba.z = __saturatef(fabs(rgba.z));
    rgba.w = __saturatef(fabs(rgba.w));
    return (uint(rgba.w * 255.0f) << 24) | (uint(rgba.z * 255.0f) << 16) | (uint(rgba.y * 255.0f) << 8) | uint(rgba.x * 255.0f);
}
__device__ float4 rgbaIntToFloat(uint c)
{
    float4 rgba;
    rgba.x = (c & 0xff) * 0.003921568627f;       //  /255.0f;
    rgba.y = ((c>>8) & 0xff) * 0.003921568627f;  //  /255.0f;
    rgba.z = ((c>>16) & 0xff) * 0.003921568627f; //  /255.0f;
    rgba.w = ((c>>24) & 0xff) * 0.003921568627f; //  /255.0f;
    return rgba;
}
//column pass using coalesced global memory reads
__global__ void
d_bilateral_filter(uint *od, int w, int h, int r)
{
    int x = blockIdx.x*blockDim.x + threadIdx.x;
    int y = blockIdx.y*blockDim.y + threadIdx.y;

    if (x >= w || y >= h)
    {
        return;
    }

    float sum = 0.0f;
    float factor = 0.0f;;
    uchar4 t = {0, 0, 0, 0};
	float tw = 0.f,tx = 0.f, ty = 0.f, tz = 0.f;
    uchar4 center = tex2D(rgbaTex, x, y);
	//t = center;
	int posIndex = 0;
    for (int i = -r; i <= r; i++)
    {
        for (int j = -r; j <= r; j++)
        {
			uchar4 curPix = {0, 0, 0, 0};
			UINT d = (UINT) sqrt((double)i*i + (double)j*j);
			if(d>r)
				continue;

			if(x + j<0||y + i<0||x + j>w-1||y + i>h-1)
			{
				factor = 0;
			}
			else
			{
				curPix = tex2D(rgbaTex, x + j, y + i);
				factor =space_weight[d] *     //domain factor
                     colorLenGaussian(curPix, center);             //range factor
			}

            tw += factor * (float)curPix.w;
			tx += factor * (float)curPix.x;
            ty += factor * (float)curPix.y;
			tz += factor * (float)curPix.z;
			sum += factor;
        }
    }
	t.w = (UCHAR)(tw/sum);
	t.x = (UCHAR)(tx/sum);
	t.y = (UCHAR)(ty/sum);
	t.z = (UCHAR)(tz/sum);
    od[y * w + x] = (((UINT)t.w)<<24|((UINT)t.z)<<16|((UINT)t.y)<<8|((UINT)t.x));
}

extern "C"
void updateGaussian(float sigma_color,float sigma_space, int radius)
{
	if( sigma_color <= 0 )
        sigma_color = 1;
    if( sigma_space <= 0 )
        sigma_space = 1;
	double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
    double gauss_space_coeff = -0.5/(sigma_space*sigma_space);  

	float color_gaussian[4*256];
	float space_gaussian[1024];

	for(int i=0;i<256*4;i++)
	{
		color_gaussian[i] = (float)std::exp(i*i*gauss_color_coeff);
		space_gaussian[i] = (float)std::exp(i*i*gauss_space_coeff);
		//if(i>100) color_gaussian[i] = 0.0f; //漫画效果
	}
// 	for(int i = -radius,int maxk=0;i<radius;i++)
// 		for(int j=-radius;j<radius;j++)
// 		{
// 			double r = sqrt((double)i*i + (double)j*j);
// 			 if( r > radius )
//                 continue;
// 			space_gaussian[maxk++] = (float)std::exp(r*r*gauss_space_coeff);
// 			//space_ofs[maxk++] = (int)(i*temp.step + j*4);
// 		}

    checkCudaErrors(cudaMemcpyToSymbol(color_weight, color_gaussian, sizeof(float)*(4*256)));
	checkCudaErrors(cudaMemcpyToSymbol(space_weight, space_gaussian, sizeof(float)*(1024)));
}

extern "C"
void initTexture(int width, int height, uint *hImage)
{
	// copy image data to array
	checkCudaErrors(cudaMallocPitch(&dImage, &pitch, sizeof(UINT)*width, height));
	checkCudaErrors(cudaMallocPitch(&dTemp,  &pitch, sizeof(UINT)*width, height));
	checkCudaErrors(cudaMemcpy2D(dImage, pitch, hImage, sizeof(UINT)*width,
		sizeof(UINT)*width, height, cudaMemcpyHostToDevice));
}

extern "C"
void freeTextures()
{
    checkCudaErrors(cudaFree(dImage));
    checkCudaErrors(cudaFree(dTemp));
}

// RGBA version
extern "C"
double bilateralFilterRGBA(uint *dDest,
                           int width, int height,
                           int radius, int iterations,
                           StopWatchInterface *timer)
{
    // var for kernel computation timing
    double dKernelTime;

    // Bind the array to the texture
    cudaChannelFormatDesc desc = cudaCreateChannelDesc<uchar4>();
    checkCudaErrors(cudaBindTexture2D(0, rgbaTex, dImage, desc, width, height, pitch));

    for (int i=0; i<iterations; i++)
    {
        // sync host and start kernel computation timer
        dKernelTime = 0.0;
        checkCudaErrors(cudaDeviceSynchronize());
        sdkResetTimer(&timer);

		dim3 blockSize(16, 16);
        dim3 gridSize((width + 16 - 1) / 16, (height + 16 - 1) / 16);

        d_bilateral_filter<<< gridSize, blockSize>>>(dDest, width, height, radius);

        // sync host and stop computation timer
        checkCudaErrors(cudaDeviceSynchronize());
        dKernelTime += sdkGetTimerValue(&timer);

    }

    return ((dKernelTime/1000.)/(double)iterations);
}

CUDA运行结果：

  • pace_sigma 和 color_sigma正向影响平滑力度，取值越大，平滑越明显。

    • pace_sigma = 10 ；color_sigma=10；

    • pace_sigma = 150 ；color_sigma=150；

  • 漫画效果

    • 在进行双边滤波时，发现有时会出现类似漫画的效果，物体的边缘有黑边

漫画效果原理：

这是由于当颜色差距过大时，值域核为0，

颜色和空间高斯值差距过大时，定义域核和值域核为0（定义域核一般不会为0）

修改代码实现漫画效果：

  • 计算高斯核数组时，void updateGaussian( ) 加入：if(i>100) color_gaussian[i] = 0.0f; //漫画效果
  • 计算颜色距离时，colorLenGaussian(uchar4 a, uchar4 b)，去掉sqrt
时间： 2025-01-02 03:59:07

双边滤波CUDA优化——BilateralFilter CUDA的相关文章

双边滤波算法的简易实现bilateralFilter

没怎么看过双边滤波的具体思路,动手写一写,看看能不能突破一下. 最后,感觉算法还是要分开 水平 与 垂直 方向进行分别处理,才能把速度提上去. 没耐性写下去了,发上来,给大伙做个参考好了. 先上几张效果图. 半径参数为10. 见图,磨皮降噪效果还不错. 具体代码如下: void bilateralFilter(unsigned char* pSrc, unsigned char* pDest, int width, int height, int radius) { float delta =

OpenCL双边滤波实现美颜功能

OpenCL是一个并行异构计算的框架,包括intel,AMD,英伟达等等许多厂家都有对它的支持,不过英伟达只到1.2版本,主要发展自己的CUDA去了.虽然没有用过CUDA,但个人感觉CUDA比OpenCL更好一点,但OpenCL支持面更管,CPU,GPU,DSP,FPGA等多种芯片都能支持OpenCL.OpenCL与D3D中的像素着色器非常相似. 1.双边滤波原理 双边滤波器的原理参考女神Rachel-Zhang的博客 双边滤波器的原理及实现. 引自Rachel-Zhang的博客,原理如下: 双

图像平滑技术之盒滤波、均值滤波、中值滤波、高斯滤波、双边滤波的原理概要及OpenCV代码实现

图像平滑是指直接对源图像的每个像素数据做邻域运算以达到平滑图像的目的.实质上主要就是通达卷积核算子实现的,卷积核算子的相关知识大家可以参考我写的博文http://blog.csdn.net/wenhao_ir/article/details/51691410 图像平滑也称为模糊或滤波,是图像处理中常用的技术之一,进行平滑处理时需要用到滤波器核(其实就是卷积核算子),根据滤波器核函数来实现不同的滤波技术.下面介绍几种 常用的图像平滑方法的大概原理及OpenCV下的实现代码. 一.盒滤波(均值滤波)

灰度图像--图像增强 双边滤波 Bilateral Filtering

学习DIP第31天 转载请标明本文出处:http://blog.csdn.net/tonyshengtan,欢迎大家转载,发现博客被某些论坛转载后,图像无法正常显示,无法正常表达本人观点,对此表示很不满意.有些网站转载了我的博文,很开心的是自己写的东西被更多人看到了,但不开心的是这段话被去掉了,也没标明转载来源,虽然这并没有版权保护,但感觉还是不太好,出于尊重文章作者的劳动,转载请标明出处!!!! 文章代码已托管,欢迎共同开发:https://github.com/Tony-Tan/DIPpro

16.双边滤波

//双边滤波 int g_nBilateralValue; void on_ChangeBilateralValue(int,void*) { if(g_nBilateralValue>0) { bilateralFilter(g_srcImage_Contrast,g_dstImage_Contrast,g_nBilateralValue,g_nBilateralValue*2,g_nBilateralValue/2); imshow("双边滤波",g_dstImage_Con

数学之路-python计算实战(18)-机器视觉-滤波去噪(双边滤波与高斯滤波 )

高斯滤波就是对整幅图像进行加权平均的过程,每一个像素点的值,都由其本身和邻域内的其他像素值经过加权平均后得到.高斯滤波的具体操作是:用一个模板(或称卷积.掩模)扫描图像中的每一个像素,用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值. #滤波去噪 lbimg=cv2.GaussianBlur(newimg,(3,3),1.8) cv2.imshow('src',newimg) cv2.imshow('dst',lbimg) cv2.waitKey() cv2.destroyAllW

【OPENCV入门之六】非线性滤波（中值滤波、双边滤波）

参考网站: http://blog.csdn.net/poem_qianmo/article/details/23184547 在很多情况下,比如在噪声是散粒噪声而不是高斯噪声时(图像偶尔会出现很大的值的时候),在这种情况下,用高斯滤波器对图像进行模糊的话,噪声是不会被去除的,它们只是转换为更为柔和但仍然可见的散粒.而用非线性滤波会更好些. 1.中值滤波(Median filter)--medianBlur函数 该方法在去除脉冲噪声.斑点噪声(speckle noise).椒盐噪声(salt-a

OpenCv高斯，中值，均值，双边滤波

#include "cv.h" #include "highgui.h" #include <iostream> using namespace std; using namespace cv; int main(int argc, char* argv[]) { Mat src = imread("misaka.jpg"); Mat dst; //参数是按顺序写的 //高斯滤波 //src:输入图像 //dst:输出图像 //Siz

【OpenCV】邻域滤波：方框、高斯、中值、双边滤波

原文:http://blog.csdn.net/xiaowei_cqu/article/details/7785365 邻域滤波(卷积) 邻域算子值利用给定像素周围像素的值决定此像素的最终输出.如图左边图像与中间图像卷积禅城右边图像.目标图像中绿色的像素由原图像中蓝色标记的像素计算得到. 通用线性邻域滤波是一种常用的邻域算子,输入像素加权得到输出像素: 其中权重核   为“滤波系数”.上面的式子可以简记为: [方框滤波] 最简单的线性滤波是移动平均或方框滤波,用 窗口中的像素值平均后输出,核函数