C语言随机半径中值滤波(扩展及百分比滤波器)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参数)。

         C语言 1         
            C语言 2

                                      图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所显示。如前所述,加法、减法以及计算直方图的中值的耗时且是有凭让图像各好的乘除,而于滤波半径无关。

                                   
  C语言 3

                图 2  算法的一定量步执行
 (a)右侧的列直方图的创新通过减最上部和加以最下像从信息获得

                     
                                                    (b)
核直方图的换代通过减最左边和加以最右面列直方图信息得到

   上述的实际效果就是核直方图为右侧走,而排直方图于下活动。在计算着,每个像素只待访问同次于,并且让上加至一个直方图数据遭到。这样,最后一步就是是计算中价值了,如代码段1所显示,这吗是一个O(1)操作。

   
综上所述,所有的单像素操作(包括更新列以及核直方图、计算中值)都是
O(1)操作。现在,我们根本来说说初始化操作,即通过积累前r行的数来算列直方图以及过去r列直方图数据测算第一单如素点的核直方图。这个历程是只O(R)操作。此外,从一行移动至下一行也占有了另外一个O(R)操作(表示未晓得)。然而,既然这初始化操作对每行只实行同样不良,相对于其它计量的话,这个好忽略不计。

                  C语言 4

    针对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的内联函数会不见面也闹立问题。

   
关于算法的耗时情况,原文为起了一个图纸:

              C语言 5

         我论机上安装的凡PS
CS4,经过测试,其中间值算法的耗时吧是乘用户输入的半径的加码而成线性增加的,但是在稍半径的时候还是比较快之。

 
   对于小半径,我之建议是采取黄算法的模式+多层直方图的法来促成,速度会比较本文要重复快来。从理论及分析,而只有当半径过7时,可应用本文效果达到O(1)复杂度。

 
 C语言 6  C语言 7  C语言 8

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

 
 C语言 9  C语言 10  C语言 11

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

     
以一副1024*768底24个真正彩色图像也例,进行速度结果反映:

时间比较(ms)
      半径         本文算法        黄算法+多层直方图
      2            506         259
      4            512         323
      6             478         377
      8            469         452
     10            479          515
     20            467         1004
     50            483         2333
     100            525         4947

 

 

 

      

 

 

 

   

 

   
同样,提供个编译好之文书为闹趣味研究该算法的心上人看效果:

 
  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

      C语言 12
     

     
由于_mm_add_epi16凡是本着16各之数据类型进行的拍卖,所以备受值得半径一般要求不盖128,否则数起数溢出等不当,工程被这样深之半径已经足足应付大部分场所的。

 

 

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