采用灰度分组统计方法实现图像中值滤波

中值滤波是图像处理中常用的一种噪声滤波方法。传统的图像中值滤波代码采用排序方法实现，处理速度主要取决于排序算法，但无论什么排序算法，总离不开大量的元素比较、交换或移动，而这些恰好是当前计算机处理的“弱项”（有经验的程序员都知道，计算机数据处理中，比较、转移、交换和频繁的数据移动比直接的算术运算和逻辑运算耗时多了），再加上没有一种好的排序算法能同时适应不同滤波半径的数据排序速度，所以在传统中值滤波实现代码中多使用选择排序、冒泡排序或者直接排序等简单排序算法，高级点的如快速排序法用在中值滤波代码中往往会使处理速度更慢。对于半径为1的中值滤波倒是有一种较好的排序算法，我在《Delphi图像处理 -- 图像的中值滤波》一文中用汇编实现过，处理速度还是较快的。

既然排序过程是图像中值滤波处理的瓶颈，能不能抛开它，用其它手段实现呢？这就是本文要探讨的问题。有朋友可能会有疑问，不排序怎么获取中间值呢，是否采用网上有些文章介绍的近似值来代替？不，本文介绍的方法决不是近似中间值，而是的的确确的“精确”中间值。

我是自学统计大专毕业，图像中值滤波中的中间值。在统计学中叫做中位数，是平均数指标的一种。平均数指标按统计的复杂程度可分为简单平均数和加权平均数，所谓简单平均数就是对统计总体的每个个体进行累计计算后求得；而加权平均数则是先对统计总体的所有个体进行分组，然后以各个组的大小作为权数后进行累计计算所得。中位数既然是平均数指标的一种，当然也存在两种计算方法。加权中位数和加权算术平均数的区别在于后者是对各个分组用权数相乘进行累积后除以总体个数，而前者只需要对各组权数进行累积后除以2，累积权数大于或等于这个数的组的值就是中位数。传统中值滤波实现代码中采用的排序手段实质就是简单中位平均数统计方法，而本文要介绍的是采用分组加权的方法来获取中位平均数，即中值滤波的中间值。

采用分组加权统计方法的前提条件是对无限统计总体进行有限的分组，例如，人口年龄统计，就是首先确定各个年龄段，这个年龄段可以是一岁、五岁、十岁等。而图像像素的R、G、B值在0 -- 255之间，正好符合有限分组的前提条件。说到这里，很多人可能已经明白我要表达的意思了：

1、按R、G、B分别定义一个256大小的灰度统计数据；

2、将图像像素总体的每个个体像素的R、G、B值进行归类；

3、确定累积中间值权数的大小；

4、从数组左端或者右端开始权数累计，累积到大于或等于中间值权数的那个灰度组就是中间值！

从前面几个步骤不难看出，前3个步骤就是图像处理中灰度统计的分类方法，第4个累积步骤与灰度统计累积计算有2个不同点，一是灰度统计累积的是权数*灰度，而这里只累积灰度；二是灰度统计累积需要全部完成，而中值累积只要达到中间值权数那个组就可以终止。在这4个步骤中，前3个步骤是相当快的（实际上只是第二个步骤），因为其中既无乘除运算，也没有比较转移，更没有元素交换或移动等动作。而制约数据处理速度的瓶颈就在第4步，因为对每个像素都必须分别按R、G、B通道对一共768个元素大小的数据进行累积，哪怕是每个通道除了2个加法运算和唯一的一次比较判断外，没有其它运算，也不需要每次都必须累积到位，但仍然是比较耗时的，不过本文在代码实现过程中，尽可能地作了一些弥补。最后结果同排序方法比起来，可算是相当快捷了，而且滤波半径越大，差距越明显，几十倍的差距绝不是天方夜谭！就算是同我前面所说的半径为1的改进排序算法比较来，大多数情况下，也略有胜出。

下面是用BCB2007写的全部实现代码（头文件略）：

//---------------------------------------------------------------------------
unsigned GetMedianGray(const BitmapData *data)
{
	unsigned buffer[256];

	memset(buffer, 0, 256 * sizeof(unsigned));
	PRGBTRIPLE p = (PRGBTriple)data->Scan0;
	unsigned width = data->Width >> 3;
	unsigned height = data->Height >> 3;
	int offset = (data->Stride << 3) - (width << 3) * sizeof(RGBTRIPLE);
	for (unsigned y = 0; y < height; y ++, (char*)p += offset)
	{
		for (unsigned x = 0; x < width; x ++, p += 8)
		{
			buffer[p->rgbtBlue] ++;
			buffer[p->rgbtGreen] ++;
			buffer[p->rgbtRed] ++;
		}
	}
	unsigned median = (height * width * 3) >> 1;
	unsigned v = 0;
	unsigned *pi = buffer;
	while (v < median)
		v += *pi ++;
	return (pi - buffer - 1);
}
//---------------------------------------------------------------------------
void GetExpendData(const BitmapData *data, unsigned radius, BitmapData *exp)
{
	exp->Width = data->Width + (radius << 1);
	exp->Height = data->Height + (radius << 1);
	exp->Stride = (exp->Width * sizeof(RGBTRIPLE) + 3) & -4;
	exp->Scan0 = (void*)new char[exp->Stride * exp->Height];
	PRGBTRIPLE ps = (PRGBTriple)data->Scan0;
	PRGBTRIPLE pd = (PRGBTriple)((char*)exp->Scan0 + radius * exp->Stride);
	PRGBTRIPLE pt = pd;
	unsigned width = data->Width;
	unsigned height = data->Height;
	int dstOffset = exp->Stride - exp->Width * sizeof(RGBTRIPLE);
	int srcOffset = data->Stride - (width - 1) * sizeof(RGBTRIPLE);
	unsigned x, y;
	for (y = 0; y < height; y ++, (char*)pd += dstOffset, (char*)ps += srcOffset)
	{
		for (x = 0; x < radius; *pd ++ = *ps, x ++);
		for (x = 0; x < width; *pd ++ = *ps ++, x ++);
		for (x = 0, ps --; x < radius; *pd ++ = *ps, x ++);
	}
	PRGBTRIPLE pb = (PRGBTriple)((char*)pd - exp->Stride);
	PRGBTRIPLE pd0 = (PRGBTriple)exp->Scan0;
	for (y = 0; y < radius; y ++, (char*)pd += dstOffset, (char*)pd0 += dstOffset)
	{
		PRGBTRIPLE p0 = pt;
		PRGBTRIPLE p1 = pb;
		for (x = 0; x < exp->Width; *pd0 ++ = *p0 ++, *pd ++ = *p1 ++, x ++);
    }
}
//---------------------------------------------------------------------------
void ImageMedianValue(BitmapData *data, unsigned radius)
{
	unsigned buffer[768];
	unsigned *pb = buffer;
	unsigned *pg = pb + 256;
	unsigned *pr = pb + 512;
	unsigned *pbOff, *pgOff, *prOff;
	int delta;
	// 获取图像总体1/8样本的灰度中间值
	if (GetMedianGray(data) < 128)
	{
		// 如果总体灰度中间值小于128，从缓冲区左边开始累计
		pbOff = pb - 1;
		pgOff = pg - 1;
		prOff = pr - 1;
		delta = 1;
	}
	else
	{
		// 否则，从缓冲区右边开始累计
		pbOff = pb + 256;
		pgOff = pg + 256;
		prOff = pr + 256;
		delta = -1;
	}

	if (radius == 0) radius = 1;
	unsigned size = (radius << 1) + 1;			// 中值滤波窗口大小
	unsigned median = (size * size + 1) >> 1;	// 缓冲区累计中值权数
	// 获取按半径扩展边框的备份图
	BitmapData exp;
	GetExpendData(data, radius, &exp);

	int stride = exp.Stride;
	int rowOffset = stride - size * sizeof(RGBTRIPLE);
	int mOffset = stride * size - sizeof(RGBTRIPLE);

	PRGBTRIPLE ps = (PRGBTriple)exp.Scan0;
	PRGBTRIPLE pd = (PRGBTriple)data->Scan0;
	unsigned width = data->Width;
	unsigned height = data->Height;
	int dstOffset = data->Stride - width * sizeof(RGBTRIPLE);
	int srcOffset = exp.Stride - (width - 1) * sizeof(RGBTRIPLE);
	unsigned i, j, v;
	unsigned *pv;

	for (unsigned y = 0; y < height; y ++, (char*)ps += srcOffset, (char*)pd += dstOffset)
	{
		// 对每行首像素按中值滤波窗口大小作一次灰度统计
		memset(buffer, 0, 768 * sizeof(unsigned));
		PRGBTRIPLE p = ps;
		for (i = 0; i < size; i ++, (char*)p += rowOffset)
		{
			for (j = 0; j < size; j ++, p ++)
			{
				pb[p->rgbtBlue] ++;		// 蓝色灰度统计
				pg[p->rgbtGreen] ++;	// 绿色灰度统计
				pr[p->rgbtRed] ++;    	// 红色灰度统计
            }
		}
		unsigned x = 0;
		while (TRUE)
		{
			// 按R、G、B灰度统计缓冲区分别进行累计，
			// 找出大于或等于median的灰度组作为中值
			for (v = 0, pv = pbOff; v < median; pv += delta, v += *pv);
			pd->rgbtBlue = pv - pb;
			for (v = 0, pv = pgOff; v < median; pv += delta, v += *pv);
			pd->rgbtGreen = pv - pg;
			for (v = 0, pv = prOff; v < median; pv += delta, v += *pv);
			pd->rgbtRed = pv - pr;
			pd ++;
			if (++ x == width) break;
			// 从灰度统计缓冲区中剔除上个像素滤波窗口最左边一列，
			// 将当前像素滤波窗口最右边一列追加到灰度统计缓冲区
			p = ps + size;
			for (i = 0; i < size; i ++, (char*)ps += stride, (char*)p += stride)
			{
				pb[ps->rgbtBlue] --;
				pg[ps->rgbtGreen] --;
				pr[ps->rgbtRed] --;
				pb[p->rgbtBlue] ++;
				pg[p->rgbtGreen] ++;
				pr[p->rgbtRed] ++;
			}
			(char*)ps -= mOffset;
		};
	}
	delete[] exp.Scan0;
}
//---------------------------------------------------------------------------

实现代码中，在前面所说的4个步骤基础上作了3点完善：

1、并非对图像的所有像素都按4个步骤进行。除了每行行首像素作了完整的分类统计外，其它像素只是从统计数组中剔除前一像素临近值的最左边一列和追加当前像素临近值的最右边一列，这无疑加快了分类统计速度，而且滤波半径越大，效果越明显。

2、在对图像进行滤波处理前，对整个图像像素总体的1/8样本求了一次中间值权数，其作用是确定像素分类后的累积方向，如果这个总体中间值权数小于128，从灰度统计数据左边（低端）开始累积，否则则从灰度统计数据右边（高端）开始累积。通过这个总体中间值权数就可以大致确定该图像处理速度的快慢，如果这个值接近两端，处理速度相对就快些，反之如靠近128附近，处理速度则相应会慢些，同样大小的图片，因这个原因，可能造成成倍的差距，但是最坏的情况下，对灰度统计数组累积时的平均值也不会超过一半。

3、制作滤波备份源图时按滤波半径对图像边缘作了扩充，这样有利于对图像边缘进行滤波处理。有些人在写中值滤波代码时往往对图像边缘略过不作处理，这在小半径滤波范围内还无所谓，但是对于较大半径的滤波处理后，由于未进行滤波处理的边缘较宽。会使得图像很难看。

正式由于第2点中所说的差距原因，本文也不好给出具体的测试数据，但有一点是肯定的：同等条件下（语言、程序员水平及编译器优化程度等），比传统排序法实现的中值滤波要快不少，滤波半径越大，差距越大！如果要求更快的处理速度，可用汇编对核心代码进行优化，我自己使用的就是如此。

最后还是贴个中值滤波测试代码和效果图：前面是源图，后面是50半径中值滤波图：

void __fastcall TForm1::Button1Click(TObject *Sender)
{
	Bitmap *bmp = new Bitmap(WideString("d:\\source.bmp"));
	BitmapData data, exp;
	Gdiplus::Rect r(0, 0, bmp->GetWidth(), bmp->GetHeight());
	bmp->LockBits(&r, ImageLockModeRead | ImageLockModeWrite, PixelFormat24bppRGB, &data);
	UINT time1 = GetTickCount();
	ImageMedianValue(&data, 50);
	Caption = GetTickCount() - time1;
	bmp->UnlockBits(&data);
	Gdiplus::Graphics *g = new Gdiplus::Graphics(Canvas->Handle);
	g->DrawImage(bmp, 0, 0);
	delete g;
	delete bmp;
}

尽管我十分努力，但水平有限，错误在所难免，欢迎指正和指导。邮箱地址：

[e-mail]