【图像处理】快速均值滤波

​忙里偷闲，整理了一下之前看的材料和自己实现的代码，这次说说均值滤波的事儿。

之前写过关于中值滤波的函数，传送门<<请戳我。

---

均值滤波，说简单点无非是一个盒子滤波器，从左上到右下划过整个图像，看似简单，可是当仔细研究，会发现，这一机械式的循环过程中包含了大量的重复计算，，，对于我这种重度强迫症患者，，，简直不能忍，，，，

---

来看一个快速的均值滤波算法思路：

因为一个数据集合的均值仅与该集合中的元素值的累加和有关,与分布形式无关,所以均值滤波采用邻域的累加值来维持交集内的像素的累加和。即:

$$average = {{sum}\over {N \times M}}$$

在同一行上,当从一个象素 $g(x,y)$ 移动到下一个象素 $g(x + 1,y)$ 时,累加和的变化为:

$$sum = sum - \sum\limits_{j=-M/2}^{M/2} g(x-\frac{N}{2},y+j) + \sum\limits_{j=-M/2}^{M/2} g(x+1+\frac{N}{2},y+j)$$

其中：
$\sum\limits_{j=-M/2}^{M/2} g(x-\frac{N}{2},y+j)$是去掉的左边列的像素的累加和；
$\sum\limits_{j=-M/2}^{M/2} g(x+1+\frac{N}{2},y+j)$是新添的右边列的象素的累加和；
在考虑换行时,相邻行之间象素的邻域也有大量的重叠区域,如何记录从一行换到下一行的交集信息呢?分析可知,实际上改变了邻域中按列的累加和，因此,还需要一个数据结构记录邻域中每列的累加和，水平滑动时去左添右，垂直换行时去上添下。

请看图加深理解

好了，有这些就够了，下面上代码

/**
 * @brief Fast mean filter
 * @param src
 * @param width
 * @param height
 * @param maskw
 * @param maskh
 */
void meanBlur(unsigned char* src, int width, int height, int maskw, int maskh)
{
    int i, j, l, temp;
    int width1 = width + maskw - 1;
    int height1 = height + maskh - 1;
    int *sumCol = NULL;
    int *psum, *ps;
    int masksz = maskw * maskh;
    unsigned char *pCurOrg, *pCurObj, *pu, *pd;
    unsigned char *src1 = NULL;
    sumCol = Malloc(int, width1);
    memset(sumCol, 0, width1 * sizeof(int));
    src1 = Malloc(unsigned char, width1 * height1);
    //获得带边界扩展的图像副本
    boundaryExtention(src, src1, width, height, maskw, maskh, width1, height1);
    //清空原图像数据
    memset(src, 0, width * height * sizeof(unsigned char));
    //构造首个横向滑窗列
    //sumCol[i] = sum(和滑窗等高的每一列像素灰度)
    for(i = 0, psum = sumCol; i < width1; i++, psum++)
        for(j = 0, pCurOrg = src1 + i; j < maskh; j++, pCurOrg += width1)
            *psum += *pCurOrg;
    for(i = 0, pCurObj = src; i < height; i++)
    {
        psum = sumCol;
        temp = 0;
        temp = sumVect(sumCol, maskw);
        pCurObj = src + i * width;
        *pCurObj++ = temp / masksz;
        for(j = 1; j < width; j++, psum++)
        {
            temp = temp - *psum + *(psum + maskw);
            *pCurObj++ = temp / masksz;
        }
        for(pu = src1 + i * width1, pd = src1 + (i + maskh) * width1, ps = sumCol, l = 0; l < width1; l++, ps++, pd++, pu++)
            *ps = *ps + *pd - *pu;
    }
    free(sumCol);
    free(src1);
}

其中，boundaryExtention()为边界扩展函数，采用复制最近元素的方式进行扩展；sumVect()是求指定长度向量的和，再次都贴出来吧。

/**
* @brief Generate a boundary-extended image to mask by boundary-replicate;
* @param src
* @param dst
* @param width
* @param height
* @param nwidth
* @param nheight
* @param new_w
* @param new_h
*/
void boundaryExtention(unsigned char* src, unsigned char *dst, int width, int height, int nwidth, int nheight, int new_w, int new_h )
{
   int i;
   int half_nw, half_nh;
   unsigned char *pCurOrg, *pCurObj;
   if(nwidth / 2 != (nwidth - 1) / 2)
   {
       printf("The mask window width is not odd !\n");
       return;
   }
   if(nheight / 2 != (nheight - 1) / 2)
   {
       printf("The mask window height is not odd !\n");
       return;
   }
   half_nw = nwidth / 2;
   half_nh = nheight / 2;
   if(new_w != width + nwidth - 1 || new_h != height + nheight - 1)
   {
       printf("The extended image space caculated is wrong!");
       return;
   }
   memset(dst, 0, new_w * new_h * sizeof(unsigned char));

   //水平对拷
   for(i = 0, pCurOrg = src, pCurObj = dst + half_nw; i < half_nh; i++)
   {
       memcpy(pCurObj, pCurOrg, width);
       pCurObj = pCurObj + new_w;
   }

   for(i = 0, pCurOrg = src, pCurObj = dst + half_nh * new_w + half_nw; i < height; i++)
   {
       memcpy(pCurObj, pCurOrg, width);
       pCurOrg = pCurOrg + width;
       pCurObj = pCurObj + new_w;
   }
   for(i = 0, pCurOrg = src + (height - 1) * width, pCurObj = dst + (height + half_nh) * new_w + half_nw; i < half_nh; i++)
   {
       memcpy(pCurObj, pCurOrg, width);
       pCurObj = pCurObj + new_w;
   }

   //垂直对拷
   for(i = 0, pCurOrg = dst + half_nw, pCurObj = dst; i < new_h; i++)
   {
       memset(pCurObj, *pCurOrg, half_nw);
       pCurOrg += new_w;
       pCurObj += new_w;
   }
   for(i = 0, pCurOrg = dst + width + half_nw-1 , pCurObj = dst + width + half_nw; i < new_h; i++)
   {
       memset(pCurObj, *pCurOrg, half_nw);
       pCurOrg += new_w;
       pCurObj += new_w;
   }
   return;
}

/**
* @brief add the vect
* @param vect
* @param length
* @return the sum
*/
int sumVect(int *vect, int length)
{
   int result=0;
   int *p,*pend;
   for(p=vect,pend=vect+length;p<pend;p++)
       result += *p;
return result;
}

经测试，该均值滤波与OpenCV和的结果都是一致的，对于512*512的lena灰度图，耗时3～4ms，，，OpenCV 开挂用了Intel 的IPP库，高富帅，比不起，，，