中值滤波原理及其C++实现与CUDA优化

对于数字图像的去噪，前边我们讲了均值滤波算法与高斯滤波算法，此外很常见的还有中值滤波算法，这些滤波算法都属于空间滤波，即对于每一个像素点，都选取其周围矩形区域中的像素点来计算滤波值。最近在项目中要使用到中值滤波，发现如果调用Opencv的medianBlur函数来实现中值滤波，窗口为3*3或者5*5时耗时为几毫秒，当窗口达到7*7或者9*9以上，耗时将增加至几十毫秒，这很影响实时性，所以自己基于C++、Opencv和CUDA实现了一个中值滤波函数，发现当窗口达到7*7以上时，比Opencv快了不少，这还是很有成就感的。所以本文的主要内容是讲中值滤波的原理、实现和优化。

1. 中值滤波的适用场景。中值滤波对盐椒噪声的去噪效果相当好，其效果比均值滤波、高斯滤波要好得多。这里的盐椒噪声，我们也可以理解为一些随机分布的、比较孤立的噪点，如下图所示。在下面的内容中，我们将使用不同的滤波算法分别对这张布满盐椒噪声的图像进行去噪，并比较去噪结果。

噪声图

2. 中值滤波原理。其原理非常简单，从字面就可以理解，所谓中值，就是对邻域矩形窗口所有点的像素值进行排序（从大到小或从小到大都ok），然后取排在最中间的像素值作为当前待滤波点的滤波值。

假设取3*3窗口进行中值滤波，如下图所示，P4为当前待滤波点，取其周围矩形区域中9个点的像素值进行排序，假设排序结果为：P0<P5<P7<P2<P3<P1<P8<P4<P6，那么则取P3的像素值作为P4的滤波值。

3. 查询中值的算法选择。查找中值的效率，是整个中值滤波算法效率的关键。所以为了尽可能地减少计算耗时，我们必须使用一种高效的查找中值算法。

(1) 冒泡排序。对矩形窗口内所有点的像素值进行从小到大或从大到小的冒泡排序，然后取排在最中间的像素值作为滤波值。这是最直接的中值查找算法，也是最耗时的算法。

冒泡排序的原理很简单，就是一轮一轮地逐个比较，假设有9个像素值P0~P8，现使用冒泡排序对其进行从大到小地排序，排序过程分为8个步骤：

步骤一：分别依次比较P0与Pi（1 ≤ i ≤ 8）的大小。

如果任意Pi大于P0，都交换Pi与P0的值，假设P3>P0，那么则交换P3与P0的值：

步骤二：分别依次比较P1与Pi（2 ≤ i ≤ 8）的大小。如果任意Pi大于P1，都交换Pi与P1的值。

步骤三：分别依次比较P2与Pi（3 ≤ i ≤ 8）的大小。如果任意Pi大于P2，都交换Pi与P2的值。

.

.

.

步骤八：比较P7与P8的大小。如果P8大于P7，则交换P8与P7的值。

(2) 快速排序。

快速排序的核心思想：首先把所有数据看成一类，选择一个界线值P0（通常把待排序数据的第一个值P0作为界线值），然后以P0为分界线，把其余大于P0的数据放到P0的左侧，小于等于P0的数据放到P0的右侧。此时P0左、右侧的数据分别组成一个新的类，然后继续对新类进行以上同样的操作，一直到所有类的数据个数都为1为止。

所以快速排序是一个裂变的过程：

(3) 分行求中值。针对于中值滤波算法，我们的目的是求取矩形窗口中所有像素值的中值，实际上只要求得中值就ok了，并不需要像上面一样对窗口中所有像素值进行排序。假设取3*3窗口，窗口内的所有像素值为P0~P8，P0~P2为窗口的一行，P3~P5为窗口的一行，P6~P8为窗口的一行，那么我们可以按照以下步骤求取3*3窗口的中值：

步骤一：对P0~P2进行排序，得到这三个数中的最小值、中值、最大值分别记为min0、mid0、max0。

步骤二：对P3~P5进行排序，得到这三个数中的最小值、中值、最大值分别记为min1、mid1、max1。

步骤三：对P6~P8进行排序，得到这三个数中的最小值、中值、最大值分别记为min2、mid2、max2。

步骤四：求min0、min1、min2的最大值max(min)。

步骤五：求mid0、mid1、mid2的中值mid(mid)。

步骤六：求max0、max1、max2的最小值min(max)。

步骤七：求max(min)、mid(mid)、min(max)的中值，即为P0~P8的中值。

这种算法虽然减少了不少计算量，但是当窗口尺寸增加之后，多次求局部最小值、中值、最大值也是比较费力的，整体耗时也就上来了。

(4) 使用窗口的灰度直方图来求取中值。

主要分为两个步骤：

步骤一：统计窗口内的灰度直方图。

步骤二：使用上一步得到的灰度直方图，从第0级灰度开始计算累加直方图。当某一级灰度i的累加直方图值达到窗口总像素的一半时，则判定i就是窗口内所有像素值的中值。

图像的灰度级是确定的，比如8位图像，其灰度级（像素值）为0~255，那么首先统计窗口内0~255像素值各自的点数，即为统计直方图，然后从0开始往后计算各个灰度级的累加直方图值，累加直方图值即对于每一灰度级，窗口内像素值小于等于该灰度级的点数。比如当计算到灰度级50时，判定发现灰度级50的的累加直方图达到窗口总像素的一半，则停止计算，取50作为中值。

综合考虑，使用直方图获取中值，计算复杂度相对更低，因此我们使用该方法来实现中值滤波。

4. 代码实现。

首先是CUDA的核函数代码：

#define WIN_SIZE 3
#define WIN_SIZE_2 (WIN_SIZE>>1)
#define D_SIZE (WIN_SIZE*WIN_SIZE)
#define LEVEL 256   //0~255


//定义纹理内存，用于快速访问源图像
texture<uchar, cudaTextureType2D, cudaReadModeElementType>  tex_src;


//源图像使用纹理内存绑定访问
__global__ void mediablur_ker(uchar *out, int row, int col)
{
  int x = blockIdx.x * blockDim.x + threadIdx.x;    //col
  int y = blockIdx.y * blockDim.y + threadIdx.y;    //row


  if (x < col && y < row)
  {
    uchar hist[LEVEL] = {0};


    for (int i = 0; i < WIN_SIZE; i++)   //统计窗口内的灰度直方图
    {
      for (int j = 0; j < WIN_SIZE; j++)
      {
        uchar idx = tex2D(tex_src, (x+j), (y+i));
        hist[idx]++;
      }
    }


    uchar mid = 0;
    uchar n = D_SIZE / 2 + 1;
    for (int i = 0; i < LEVEL; i++)   //计算累加直方图
    {
      mid += hist[i];
      if (mid >= n)  //当某一级灰度i的累加直方图值达到窗口总像素的一半时，则判定i就是窗口内所有像素值的中值
      {
        out[y*col + x] = (uchar)i;
        break;
      }
    }
  }
}

其次是实现函数代码：

void mediablur_fiter_cuda(Mat src, Mat &dst)
{
  Mat src_board;
  //扩充边缘
  copyMakeBorder(src, src_board, WIN_SIZE_2, WIN_SIZE_2, WIN_SIZE_2, WIN_SIZE_2, BORDER_REFLECT);   //扩充边缘


  const int row0 = src.rows;
  const int col0 = src.cols;
  const int img_size0 = row0*col0;


  const int row = src_board.rows;
  const int col = src_board.cols;
  const int img_size = row*col;


  //////////////////////////////////////////////////////////////////////////////////////////////////////////////
  uchar *dst_cuda;
  cudaMalloc((void**)&dst_cuda, img_size0);
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////


  cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<uchar>();//声明数据类型
  cudaArray *cuArray_src;
  cudaMallocArray(&cuArray_src, &channelDesc, col, row);  //分配大小为col*row的CUDA数组
  tex_src.addressMode[0] = cudaAddressModeWrap;//寻址方式
  tex_src.addressMode[1] = cudaAddressModeWrap;//寻址方式 如果是三维数组则设置texRef.addressMode[2]
  tex_src.normalized = false;//是否对纹理坐标归一化
  tex_src.filterMode = cudaFilterModePoint;//纹理的滤波模式：最近点取样和线性滤波  cudaFilterModeLinear
  cudaBindTextureToArray(&tex_src, cuArray_src, &channelDesc);  //纹理绑定，CUDA数组和纹理参考的连接


  cudaMemcpyToArray(cuArray_src, 0, 0, src_board.data, img_size, cudaMemcpyHostToDevice);


  //////////////////////////////////////////////////////////////////////////////////////////////////////////////


  dim3 Block_G(16, 16);
  dim3 Grid_G((col0 + 15) / 16, (row0 + 15) / 16);
  //调用核函数
  mediablur_ker << <Grid_G, Block_G >> >(dst_cuda, row0, col0);
  
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////
  dst = Mat::zeros(src.size(), CV_8UC1);
  cudaMemcpy(dst.data, dst_cuda, img_size0, cudaMemcpyDeviceToHost);


  //cout << dst << endl;
  //////////////////////////////////////////////////////////////////////////////////////////////////////////////


  cudaFree(dst_cuda);
  cudaFreeArray(cuArray_src);
  cudaUnbindTexture(tex_src);
  
}

5. 滤波结果对比。

(1) 统一取3*3窗口，分别调用以上实现的中值滤波函数，和Opencv的高斯滤波函数、均值滤波函数对本文开头的噪声图进行滤波去噪，得到的结果如下图所示。可以看到，中值滤波对盐椒噪声的去噪效果，相比于高斯滤波和均值滤波来说，那是相当出色。

均值滤波结果

高斯滤波结果

中值滤波结果

(2) 比较以上实现的中值滤波与Opencv实现的中值滤波函数的耗时。如下表所示，可以看到，当窗口在7*7以上时，本文的优化加速效果就很明显了。