任意半径中值滤波（扩展至百分比滤波器）O(1)时间复杂度算法的原理、实现及效果。

主要参考论文：Median Filter in Constant Time.pdf

参考代码：http://files.cnblogs.com/Imageshop/CTMF.rar

中值滤波是一种经典的图像操作，特别适用于椒盐噪音的去除。同样，他也是USM锐化（表示怀疑，我记得是高斯滤波)、顺序处理、形态学操作（比如去孤点）等算法的基础。更高级别的应用包括目标分割、语音和文字识别以及医学图像处理等。

然而，过多的处理时间严重的限制住了中值滤波器的使用。由于其算法的非线性和不可分离性普通的优化技术并不合适。最原始的步骤就是获取图像像素的一个列表，然后进行排序，接着取中值。在一般的情况下，该算法的复杂度是O(r2logr)，其中r为核的半径。当像素的可能取值是个常数时，比如对于8位图像，可以使用桶排序(Bucker sort），该排序使得算法的复杂度降低为O(r2)。但是除了小半径的情况外，这样的改进任然是不可接受的。

这里插一句，从我个人的认知上说，任何基于排序的中值滤波，都是无法对大半径进行实时有效的处理的。

在<A Fast Two-Dimensional Median Filtering Algorithm>一文中，提出了基于直方图统计的快速中值滤波，其时间复杂度为O(r)，这个想法我在我的Imageshop软件里也曾经实践过（那个时候可没有看过这个论文，可见这个优化是个大众都能想到的一步）。后来发现，在Paint.net里的中值用的也是这个。具体来讲，这个算法通过计算在核内部的像素的直方图间接得到中值。当从一个像素移动到另外一个与之相邻（水平或垂直方向）的像素时，只需要更新一部分的信息，如图1所示。为了更新直方图，2r+1次加法以及2r+1次减法需要执行，而从直方图中计算中值所需要的时间是一定的，如代码段1所示。对于8位图像，直方图由256个元素组成，在平均上说，计算中值需要128次比较和127次加法。实际上，通过改变终止寻找的条件我们可以计算任何其它百分比效果（见代码段1中的Percentile参数）。

图1 黄氏算法，本图半径r=2

int StopAmount = (2 * Radius + 1) * (2*Radius +1) * Percentile;
int Sum = 0;
for (int I = 0; I <= 255; I++)
{
    Sum += HistBlue(I);
    if (Sum >= StopAmount)      // 满足停止条件
    {
        Value = I;              // 得到中值
        break; 
    }
}

代码1： 从直方图中计算中值（Percentile=0.5）或任意百分比值

　　为了使中值滤波的时间复杂性降低至线性以下，人们做出了很多努力。Gel使用了基于树的算法将复杂度降低为O(log2r)，在同一篇论文中，他们简要的说了一种复杂度为O(log r)的二维图像中值算法。我们这里提出的算法也是用直方图的计算来取代如龟速的排序。最近，Weiss使用层次直方图技术获取了O(log r)复杂度的效果，在他的方法中，虽然速度是快了，但是代码复杂难懂。本文提出一种简单而有效的算法，易于用CPU或其他硬件实现。

为更好的理解文章的算法，我们先来看看黄氏算法的不足。特别注意到该算法行与行之间没有任何信息得到保留，而每个像素的处理至少有2r+1次加法和减法的直方图计算，这就是其复杂度为O(r)的原因。凭直觉，我们猜想应该有某种方法使得对每个像素，直方图只需累加一个固定的次数，从而获得O(1)的复杂度。正如我们所看到的，通过保留行与行之间的信息，这变得可行。首先让我们来介绍下直方图的一些属性上。对于不相邻的区域A和B，有下式成立：

H(A ∪ B) = H(A) + H(B)

　　注意到两个直方图的累加是一个O(1)操作。他和直方图的元素个数有关，而直方图元素个数是由图像位深决定的。有了这一点，我们就可以开发一个新的算法。

首先，对于每一列图像，我们都为其维护一个直方图（对于8位图像，该直方图有256个元素），在整个的处理过程中，这些直方图数据都必须得到维护。每列直方图累积了2r+1个垂直方向上相邻像素的信息，初始的时候，这2r+1个像素是分别以第一行的每个像素为中心的。核的直方图通过累积2r+1个相邻的列直方图数据获取。其实，我们所做的就是将核直方图分解成他对应的列直方图的集合，在整个滤波的过程中，这些直方图数据在两个步骤内用恒定的时间保持最新。

考虑从某个像素向右移动一个像素的情况。对于当前行，核最右侧的列直方图首先需要更新，而此时该列的列直方图中的数据还是以上一行对应位置那个像素为中心计算的。因此需要减去最上一个像素对应的直方图然后加上其下面一像素的直方图信息。这样做的效果就是将列直方图数据降低一行。这一步很明显是个0(1)操作，只有一次加法和一次减法，而于半径r无关。

第二步更新核直方图，其是2r+1个列直方图之和。这是通过减去最左侧的列直方图数据，然后再加上第一步所处理的那一列的列直方图数据获得的。这一步也是个O(1)操作，如图2所示。如前所述，加法、减法以及计算直方图的中值的耗时都是一些依赖于图像位深的计算，而于滤波半径无关。

　图 2 算法的两步执行 (a)右侧的列直方图的更新通过减最上部和加最下部像素信息获得

(b) 核直方图的更新通过减最左侧和加最右侧列直方图信息获得

　 上述的实际效果就是核直方图向右移动，而列直方图向下移动。在计算中，每个像素只需访问一次，并且被添加到一个直方图数据中。这样，最后一步就是计算中值了，如代码段1所示，这也是一个O(1)操作。

综上所述，所有的单像素操作（包括更新列以及核直方图、计算中值）都是 O(1)操作。现在，我们重点来说说初始化操作，即通过累积前r行的数据来计算列直方图以及从前r列直方图数据计算第一个像素点的核直方图。这个过程是个O(R)操作。此外，从一行移动到下一行也占用了另外一个O(R)操作(表示不理解)。然而，既然这个初始化操作对每行只执行一次，相对于其他计算来说，这个可以忽略不计。

针对8位灰度图像，我们对上述算法进行一下总结。

（1）、对核最右侧的列直方图执行一次加法。

（2）、对同一列直方图执行一次减法，去除多余的像素信息。

（3）、将更新后的列直方图数据加到核直方图中，这进行了256次加法。

（4）、将无效的列直方图数据从核直方图中减去，这需要256次减法。

（5）、为找到核直方图的中值，平均需要128次比较和127次加法。

上述计算量看起来比较多。然而，大部分操作都可并行处理，从而明显的降低了处理时间。更重要的，还有很多优化措施能降低计算量，如下所示。

1、并行化

　　现代处理器提供了SIMD指令可以用来加速我们的算法。上面描述的操作大部分都是对直方图数据进行加和减操作。通过MMX,SSE2或Altivec指令可以并行处理多个直方图操作。为了在一条指令中做更多的直方图处理，直方图的只能用16位的数据类型。因此，核心的大小最大为216(对应的核的大小是（256/2-1))，对于一般的应用足够了。这个限制并不只是对我们的算法：这是一种优化的手段而已。

另外一个可以运行并行的地方就是从图像中读取数据以及将其累加到对应的直方图中。同上述交替更新列和核直方图不同的是，我们可以首先更新整行的列直方图。然后利用SIMD指令，我们可以并行的更新多列直方图数据（不同列直方图的更新之间没有任何关系，所以好并行）。然后再像平常一样更新核直方图。

这条优化手段对于有些高级语言是无法实现的。像VC6，VC.NET这类可以直接内嵌汇编的语言，虽然可以实现，也需要作者具有很好的汇编语言基础，因此，实施的难度比较大。有兴趣的读者可以参考附件中的中的SSE代码。

2、缓存优化

　　恒常时间的中值滤波算法需要在内存中为每列保持一个直方图，对于图像，这很容易就多达数百KB的大小，通常这大于今天的处理器的缓存。这导致访问内存的效率降低。一种方式就是将图像在水平方向上分成几部分分开处理。每个块的宽度大小需精心选择， 使得其对应的列直方图不大于缓存的大小而有能充分利用缓存。这种修改的不利点就是扩大的边缘效应。在实践中，这通常能导致处理时间的大幅降低。请注意，在不同的处理器上同时处理这些块是该算法的一种很简单的并行算法。

这种优化说实在我不知道如何用代码去实现。

　　3、多层直方图

　　在<A Coarse-to-Fine Algo-rithm for Fast Median Filtering of Image Data With a Huge Number of Levels>一文中显示了多尺度直方图是非常有效的优化手段。其想法是维持一个平行的较小的直方图，直方图记录了图像的高位数据。例如，对于8位图像，使用两层的直方图很常用，其中高层使用4位，而低层使用全8位数据。习惯上我们分别给他们命名为粗分和细分直方图。粗分直方图包含了16(2^4)个元素,每个元素是对应的细分直方图高位的累积和。

使用多层直方图有两个好处，第一个就是计算中值过程的加速。我们可以首先在粗分数据中需找到中值在细分数据中段的位置而不用检查整个256个位置。平均上说这只需要16次而不是128次比较和相加。第二个好处是关于直方图的相加和相减。当对应的粗分数据为0时，则可以不用计算对应的细分段。当半径 r很小时，列直方图是稀疏分布的，这个时候的分支判断是很有必要的。

以上说的很笼统。举个简单的例子吧，对于3*3的一个核，假如像素的数据值分别为： 100、120、98、77、215、50、243、199、180。对应的高位序列为：6、7、6、4、13、3、15、12、11。低位序列为：4、8、2、13、7、2、3、7、4。 那么粗分直方图的16个元素的值分别为:

Coarse[3]=1　　Coarse[4]=1　　Coarse[6]=2　　Coarse[7]=1　　Coarse[11]=1　　Coarse[12]=1　　Coarse[13]=1　　Coarse[15]=1，其他都为0；

　　中位数的累加值为3*3/2=5，对粗分直方图进行累加：Coarse[3]+Coarse[4]+Coarse[6]+Coarse[7]=5，得到中值对应的高位段是H=7，因此，中间值就应该在112-128之间，然后再从细分直方图的对应段中按照类似的方式找到低位的值L，最后得到中值H*16+L。

4、条件更新核

粗分和细分直方图的分离还有一个不明显但是很有效的优化效果。注意到算法的大部分时间都耗费更新在核直方图的时加上或减去列直方图数据，这个时间随着实时更新粗分直方图而有条件的更新细分直方图而得到降低。

记得前面说过计算中值的过程是先在粗分数据中寻找中值所在段，然后再从细分数据中找到精确值。对于核的中值，每个列直方图最多只会有2r+1次贡献，意味着只有2r+1个对应的细分段对计算结果有用。那些从来未被使用的段，其对应的细分数据将无需更新。

　 为了实现该功能，我们需要为每个开辟一个记录其最后被更新的位置的列表。当从一个像素移向下个一个像素时，我们更新列直方图以及核直方图的粗分数据。然后根据粗分数据计算出中值再细分数据中所在的段。下一步，根据这个段上次被更新的位置更新的细分直方图。如果上次更新的位置和当前列的位置相差2r+1的距离，那说明旧的位置和当前位置没有任何交叉。因此我们从0开始更新段。正是通过这种方式，我们加速了整个处理速度。

交错布置直方图数据，从而使得相邻列的直方图数据在内存也是相邻的是有好处的。因此，细分直方图数据需要按下述方式布置：段索引、列索引、最后是元素索引。这样，核的细分直方图的更新就是对一块连续的内存的累加了，具体的讲，细分直方图有类似如下的定义形式:int [,,] Fine= new int [16,Width,16],其中第一个16对应段索引，即像素值的高4位，最后一个16是元素值，即像素的低4位值。因为三维数组的访问会造成冗余的计算下标的过程，因此，为了提高速度，应该使用一维数组或者直接用指针访问内存，以一维数组为例，此时的定义应该修改为int [] Fine=new int[16*Width*16],对于第X行某个像素值Value，其对应的细分位置则为： Fine[(Value>>4) * (Width*16)+X*16+ (Value & 0xF)];

具体的可能还是看参考代码更容易理解。

参考代码中如下函数：

static inline void histogram_sub( const uint16_t x[16], uint16_t y[16] )
{
    int i;
    for ( i = 0; i < 16; ++i ) {
        y[i] -= x[i];
    }
}

第一个建议是把这个小循环手工展开。第二，我是是用C#编程实现结果的，C#没有inline，于是我把这样的代码直接展开内嵌到我的代码中，可是令人诧异的结果却是调用函数的版本速度却比内嵌的速度还要快，反汇编函数版本代码如下：

  y[0] -= x[0];
00000000  push        ebp 
00000001  mov         ebp,esp 
00000003  push        esi 
00000004  mov         esi,ecx 
00000006  mov         ecx,edx 
00000008  mov         eax,dword ptr [esi] 
0000000a  sub         dword ptr [ecx],eax 
            y[1] -= x[1];
0000000c  lea         edx,[ecx+4] 
0000000f  mov         eax,dword ptr [esi+4] 
00000012  sub         dword ptr [edx],eax 
            y[2] -= x[2];
00000014  lea         edx,[ecx+8] 
00000017  mov         eax,dword ptr [esi+8] 
0000001a  sub         dword ptr [edx],eax 

...........................其他的减法.....................................

            y[14] -= x[14];
00000074  lea         edx,[ecx+38h] 
00000077  mov         eax,dword ptr [esi+38h] 
0000007a  sub         dword ptr [edx],eax 
            y[15] -= x[15];
0000007c  add         ecx,3Ch 
0000007f  mov         edx,ecx 
00000081  mov         eax,dword ptr [esi+3Ch] 
00000084  sub         dword ptr [edx],eax 
00000086  pop         esi 
        }
00000087  pop         ebp 
00000088  ret 

关于这个问题，我的分析是，写成函数的版本，虽然多了几句压栈和出栈的语句外，CPU会充分利用寄存器来进行操作，而内嵌后，由于变量太多，CPU只能利用内存来处理这些来处理这些赋值的。而寄存器比内存的操作快多了。因此不晓得VC的内联函数会不会也有这问题。

关于算法的耗时情况，原文给出了一个图表：

我本机上安装的是PS CS4，经过测试，其中间值算法的耗时也是随着用户输入的半径的增加而成线性增加的，但是在小半径的时候还是比较快的。

对于小半径，我的建议是采用黄算法的模式+多层直方图的方式来实现，速度会比本文要更快些。从理论上分析，而只有在半径大于7时，可采用本文效果达到O(1)复杂度。

　　　　　　　　　　原始图像　　　　　　　　　　　　　　　　　　半径=5，百分比=50　　　　　　　　　　　　　　 半径=40，百分比=50

半径=5，百分比=25 　　　　　　　　　　　　 半径=5，百分比=75 半径=40，百分比=75

以一副1024*768的24位真彩色图像为例，进行速度结果汇报：

同样，提供个编译好的文件给有兴趣研究该算法的朋友看看效果：

http://files.cnblogs.com/Imageshop/MedianFilter.zip

由于上述MedianFilter是用VB6编制的，而VB6不支持指针，只能用数组来操作图像数据，在有些地方会造成效率的严重丢失，我用C#写的速度（未用多线程）是该MedianFilter的速度的三倍到4倍左右，如果是处理灰度，基本上又可一达到同样大小彩色图像速度的2到3倍。

在实际应用中，大半径的中值对于去除椒盐噪音是没有意义的，因为此时图像已经损失了太多有用的信息了。根据我的了解，大半径可以发挥用处的地方有：1、如果你的程序有和PS一样的选区技术，那么选区的平滑这个功能其实就是对选区数据进行中值处理的过程，这个当然希望之星速度和半径无关。2、一些二值图像的去噪上可以利用，比如一定半径的中值，可以去除一些孤立的块，从而消除噪音。3、在对某些图像进行艺术处理时，需要大半径的中值。

2013.10.16 补充：

最近开始学习了下C，于是用C语言实现上述过程。而在论文作者的原始代码中，用到SSE2的相关函数进行处理，即有如下的代码：

__inline void  HistgramAdd( unsigned short *x, unsigned short *y  )
{
    *(__m128i*) y = _mm_add_epi16( *(__m128i*) y, *(__m128i*) x );
    *(__m128i*) (y+8) = _mm_add_epi16( *(__m128i*) &y[8], *(__m128i*) &x[8] );
}

__inline void  HistgramSub(unsigned short *x, unsigned short *y )
{
    *(__m128i*) &y[0] = _mm_sub_epi16( *(__m128i*) &y[0], *(__m128i*) &x[0] );
    *(__m128i*) &y[8] = _mm_sub_epi16( *(__m128i*) &y[8], *(__m128i*) &x[8] );
}

这里_mm_add_epi16是一组Instructions指令，编译器在编译的时候会将其编译为SSE对应的汇编代码，而由于这些指令能实现指令级别并行，比如上述_mm_add_epi16可以在同一个指令周期对8个16位数据同时进行处理，并且HistgramAdd这些函数在程序里会大量使用到，因此程序的速度能大幅提高。

对于上述测试的1024*768的图片，用C编写DLL，C#调用，同样的机器耗时平均在160ms左右，速度较纯粹的C#的代码快了3倍多，可见SSE的强大。

[DllImport("MedianBlur.dll", CallingConvention = CallingConvention.StdCall, CharSet = CharSet.Unicode, ExactSpelling = true)]
private static extern void MedianBlurGray( byte * Src, byte * Dest, int Width, int Height ,int Stride ,int Radius,int Percent);
[DllImport("MedianBlur.dll", CallingConvention = CallingConvention.StdCall, CharSet = CharSet.Unicode, ExactSpelling = true)]
private static extern void MedianBlurRGB( byte * Src, byte * Dest,int Width, int Height ,int Stride ,int Radius,int Percent);

附件为C#调用C编写的DLL的工程代码文件：http://files.cnblogs.com/Imageshop/MedianBlurTest.rar

由于_mm_add_epi16是针对16位的数据类型进行的处理，所以中值得半径一般要求不大于128，否则数出现数据溢出等错误，工程中这么大的半径已经足够应付大部分场合的。

***************************作者： laviewpbt 时间： 2013.4.26 联系QQ: 33184777 转载请保留本行信息*************************