忙里偷闲,整理了一下之前看的材料和自己实现的代码,这次说说均值滤波的事儿。
之前写过关于中值滤波的函数,传送门<<请戳我。
均值滤波,说简单点无非是一个盒子滤波器,从左上到右下划过整个图像,看似简单,可是当仔细研究,会发现,这一机械式的循环过程中包含了大量的重复计算,,,对于我这种重度强迫症患者,,,简直不能忍,,,,
来看一个快速的均值滤波算法思路:
因为一个数据集合的均值仅与该集合中的元素值的累加和有关,与分布形式无关,所以均值滤波采用邻域的累加值来维持交集内的像素的累加和。即:
average=sumN×M
在同一行上,当从一个象素
g(x,y)
移动到下一个象素
g(x+1,y)
时,累加和的变化为:
sum=sum−∑j=−M/2M/2g(x−N2,y+j)+∑j=−M/2M/2g(x+1+N2,y+j)
其中:
∑j=−M/2M/2g(x−N2,y+j)
是去掉的左边列的像素的累加和;
∑j=−M/2M/2g(x+1+N2,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;
}