何凯明在导向滤波一文的有关材料中提供了其matlab代码,何凯明在导向滤波一文的连锁资料中提供了其matlab代码

  注意最终的IM_ClampFHtoByte函数是将括号内的值限制在0和255里边的。

 for (int Y = 0; Y < SmallH; Y++)
 {
     float PosX = 0;
     float AddX = (SmallW - 1.0f) / Width;        //  主要是为了减少下面插值向右增1的指针超过范围,但这样做其实是和精确的算法有一点点差异的
     float *LinePDA = TempA + Y * Width * Channel;    //  TempA和TempB为临时分配的大小为(SmallH * Width * Channel * sizeof(float)大小的内存
     float *LinePDB = TempB + Y * Width * Channel;
     float *LinePA = MeanA + Y * SmallW * Channel;
     float *LinePB = MeanB + Y * SmallW * Channel;
     if (Channel == 1)
     {
         for (int X = 0; X < Width; X++)
         {
             int XX = (int)PosX;
             float PartXX = PosX - XX;
             float InvertXX = 1 - PartXX;
             float *PtLeftA = LinePA + XX;
             float *PtLeftB = LinePB + XX;
             LinePDA[X] = PtLeftA[0] * InvertXX + PtLeftA[1] * PartXX;
             LinePDB[X] = PtLeftB[0] * InvertXX + PtLeftB[1] * PartXX;
             PosX += AddX;
         }
     }
   //  ...................
 }

  用_mm_loadl_epi64加载8个字节数据到内存,并用_mm_cvtepu8_epi16将其转移为16位的__m128i变量,后边再把低4位和高4位的数码分别转换成32位的,然后用_mm_add_epi32增进,注意前面我改换函数用了二种差别的主意,因为此处的数量相对都是正数,由此也是可以使用_mm_cvtepi16_epi32和_mm_unpackhi_epi16组合Zero实现。

 

     

  这么些的优化对于BoxBlur来说是重中之重的一步,不过更关键的是他得以采用在重重场面,比如图像的部分均方差计算,也足以使用类似的技能拓展加快,两幅图像大的有的平方差也是可以这么优化的,后续我会不难的谈下他在那地方加速的施用。

   
 首先,第一,步骤6中的八个采样进度不要分开写,直接写到同一个for循环内部,那样可以节省不可胜数坐标的臆想进度,第二,那里一般的上采样平常拔取双线性插值就OK了,互联网上有很多关于双线性插值的SSE优化的代码,可是那多少个基本都是本着32位的图像做的优化,搬到24位和8位中是不适用的,而我辈会在50%以上的概率中相遇24位图像,所以说啊,网络上的东西虽多,但精华太少。

     
使用SSE优化能将上述进度提速2倍以上。

           ca88亚洲城网站 1   
 ca88亚洲城网站 2

      何凯明在2015又刊出了一篇《法斯特Guided Filter》的稿子,演说了一种很实用的更快速的导向滤波流程:

     
那样改动后,所有的boxfilter均是对下取样后的数据开展拍卖,当s=4时,总括量裁减到原来的1/16,而s=5,则缩减到了原来的1/25,当然那些时候多了2个下取样和2个上取样的算法,下取样由于是压缩,总结量很小,无需关切,而上采样,总结量和原图大小有关,按照自家的估测,那些上采样的耗时可能占所有经过的一般时间左右的,是不行值得注意优化的。

                       
  ca88亚洲城网站 3

   
 由于在总结进程中真的存在有的结实大于了0和255的界定,因而只要把IM_ClampFHtoByte函数去除,对有些图像会油可是生噪点,因而,我们无法一心看重编译器的向量化优化了,那就非得团结写SIMD指令,由于SIMD自带了饱和处理的连锁函数,而上述内部的X
的for循环是很不难用SSE处理的,唯一必要小心的就是索要把LinePS对应的字节数据转换为浮点数据,那里我概括的升迁可以用如下指令将8个字节数据转换为8个浮点数:

     
可是就是是如此,由于6次统计以及中间的别样一些浮点运算,依旧给全部算法带来了很大的演算费用和内存费用,在多如牛毛场所依然不能够满意必要的,比如实时去雾等现象。在最初我的敏捷去雾已毕中,都是先使用下采样图的导向滤波结果,然后再双线性插值放大获得大图的透射率图,即使在视觉效果上能化解去雾算法的进程难点,可是若是是其他场景的导向滤波需求,依旧会合到众多欠缺的。

  往日认为那一个算法难以使用SSE优化,首要难题就在此地,可是在读书了Opencv的积分图技术时,那几个标题就缓解了,进一步分析还发现Opencv的代码写的并不健全,还有革新的空间,见上述代码,使用两遍对同一数据移位就足以拿走丰盛,由P3
P2 P1 P0变为P3+P2+P1+P0 P2+P1+P0 P1+P0 P0。我们略微分析一下。          
  

 

   
 本文纯属计流水账,未做详细分析。

  再来看看92到95行代码,那一个也很简单。

   
 在一台I5的台式机上,选取默许设置,以本人为导向图处理3000*2000的24位图像必要约55ms,假如是灰度图差不多是20ms,那个和优化后的
boxblur速度基本一致,若是开启动多线程,比如开七个线程,还是能提速25%左右,再多也无帮忙了。

 

 

       
 ca88亚洲城网站 4 
 ———–——>   ca88亚洲城网站 5 
 ———–——>  ca88亚洲城网站 6

 

  那里也是一回性处理8个像素,那里自己利用了别的一种转移技术来把8位转换为16位,可是前面的Diff因为有正有负,要转换为32位就务须利用_mm_cvtepi16_epi32来转换,无法用unpack序列组合函数来兑现,因为unpack不会活动符号位,我在此地吃了一些次亏了。

   
 本文纯属计流水账,未做详细分析。

 

 

     
有不少有情人可能不了解,假设把上边的IM_ClampFHtoByte这一个函数去掉,直接运用括号内的代码,VS的编译器可以很好的对地点代码举行向量化编译(VS编译只要你没有把代码生成–》启用增强指令集设置成无增强指令/arch:IA32,哪怕设置为未安装,都会把浮点的代码编译为SIMD相关指令的),而如若大家对分化的Channel,比如3通道4坦途在循环里展开后,很衰颓,按照大家的举行循环的争论,速度应该加紧,但实际却反倒了。所以大家需要充裕了然编译器的向量化特性,就能写成更迅捷的代码。

     
自以为当前本人优化的快慢在CPU版本中很难有人能超越了(仅仅使用CPU、不用多线程,下采样率0.2),倘使什么人有更快的算法,在第三方公证的情景下,我甘愿提供1000元奖励^_^。

 

               R = R1 * (1 – 0.4) + R2
* 0.4;

  而由第一个图到第一个图的长河大致可实用上边的代码表述:

ca88亚洲城网站 7

     
自认为方今自己优化的快慢在CPU版本中很难有人能跨越了(仅仅使用CPU、不用二十四线程,下采样率0.2),假如哪个人有更快的算法,在第三方公证的情景下,我情愿提供1000元奖励^_^。

     
我使用的一个优化措施时,先进行水平方向的上采样到一个缓冲区中(Width  *
SmallH),然后在从那几个缓冲区中沿着中度方向缓冲到(Width *
Height),如下图所示:

     
SSE图像算法优化序列五:超高速指数模糊算法的完毕和优化(10000*10000在100ms左右贯彻) 一文中,我已经说过优化后的ExpBlur比BoxBlur还要快,这么些时候我相比的BoxBlur算法是通过积分图+SSE完毕的,我在09年此外一个博客账号上早已提供过一篇这些小说彩色图像高速模糊之懒惰算法,里面也介绍了一种高效的图像模糊算法,那些算法的执行时间基本也是和半径毫不相关的。在今年的SSE优化学习之路上自家早已也考虑过将该算法使用SSE达成,但当时认为那一个算法逐像素同时逐行都是上下着重的(单纯的逐像素看重算法我在指数模糊里有关联怎么样用SSE优化),是力不从心用SSE处理的,一向没考虑,直到日前有意中人提议某个基于局地局方差的算法希望能提速,我又重新触发灵感,终于将那种算法也促成的一声令下集完成,并且测试速度比积分图快二倍,比解析opencv中Box
Filter的完毕并提议越来越加快的方案(源码共享)
这篇小说的速度快3倍,比opencv的cvSmooth函数快5倍,在一台老旧的I3台式机上拍卖3000*2000的灰度图达到了6ms的进程,本文分享该优化进程并提供灰度版本的优化代码供大家学习和议论。

     
别的一个标题,在上头的流程2的第一步中,对boxfilter的半径r也是展开了同比例的减弱的,注意到boxfilter的半径日常情况下大家都是用的平头,要是减弱后的r’也开展取整的话,举例来说,对于s
=4的图景下,半径为8、9、10、11那多样情景最后收获的导向滤波结果就全盘相同了,就如那不符合大家对算法严酷性的必要,所以我们要扶助一种浮点半径的boxfilter。

       7:   q = fupsample(q,
s)

  在彩色图像高速模糊之懒惰算法一文附带的代码中(VB6的代码)是针对24位的图像,为了啄磨方便,大家先贴出8位灰度的C++的代码:

ca88亚洲城网站 8

     
很明显,由于I的参与统计,何的做法能更大程度上维持结果和原汁原味的近乎,而自我的不二法门则会发出较大的块状相似,所以住户大神就是大神。

  接着是FillLeftAndRight_Mirror_C函数的优化,改写如下:

     
 由地方的第二个图到首个图的光景代码如下:

     
里面的浮点计算的历程的SSE代码就和平凡的函数调用没什么却别,最终的写到LinePD这一个字节数据的长河可以用_mm_storel_epi64以及有关活动搞定。

  原理基本上就是这样,这个算法占用的额外内存很明显很少,但是不支持In-Place操作。
  为了追求速度,我们把整个过程能用SSE优化的地方都用SSE优化。
  首先是第79至86行的数据,这个很容易优化:

for (int Z = -Radius; Z <= Radius; Z++)
{
    unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride;
    int BlockSize = 8, Block = Width / BlockSize;
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        int *DestP = ColValue + X + Radius;
        __m128i Sample = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(LinePS + X)));
        _mm_storeu_si128((__m128i *)DestP, _mm_add_epi32(_mm_loadu_si128((__m128i *)DestP), _mm_cvtepi16_epi32(Sample)));
        _mm_storeu_si128((__m128i *)(DestP + 4), _mm_add_epi32(_mm_loadu_si128((__m128i *)(DestP + 4)), _mm_unpackhi_epi16(Sample, _mm_setzero_si128())));
    }
    //  处理剩余不能被SSE优化的数据
}

     
那里如此做的其余一个好处是在Y循环中计算是单独的,因而都得以动用OPENMP加速。

  本文简要的记录了自身在优化导向滤波完成的经过中所适用的优化措施和一部分细节,避防时间久了后自己都不记得了,不过请不要向自己直接索取源代码。

  至于24位图像的优化,是个相比窘迫的情形,因为SSE三次性处理4个32位数,而24位各种像素唯有3个轻重,那种情景相似仍然把她恢弘到4个轻重,比如说ColValue就足以改成4坦途的,然后累积和也亟需处理成4通路的,这本来须要一定的迷你淫技,这里自己不想多谈,留点东西给自己。当然也足以设想先把24位分解到3个灰度内存,然后使用灰度的算法进行测算,后边在合成到24位,这几个解释和合成的历程也是足以行使SSE加快的,假诺你结合三十二线程,还能把3个灰度进程的处理并行化,那样除了多占用内存外,速度比至二级处理24位要快(因为一直处理算法不可以并行的,前后器重的案由)。其余就是终极在测算列累积求平均值的进度也变得越来越自然了,不会出现灰度这种要在__mm128i中间举行添加的经过,而是径直得三个SSE变量累加。

 

  自从何凯明提议导向滤波后,因为其算法的简单性和管事,该算法获得了常见的施用,以至于新版的matlab都将其看做规范自带的函数之一了,利用她可以解决的具备的保边滤波器的能缓解的标题,比如细节增强、HDR压缩、细节羽化、去雾、风格化,而且由于其保边特性,假如过多观念函数中利用高斯滤波或者均值滤波的地点用她取代,能很好解决一些强边缘的对接不自然难点,比如retinex、Highlight/shadow等使用中,由此,飞快的贯彻该算法具有很强的适用意义。

  我们定义一个Radius + Width + Radius的内存数据ColValue,用来保存列方向的累加数据,对于第一行数据,需要做特殊处理,也就是用原始的方式计算一行像素所有元素在列方向上的和,
详见78至于86行代码,当然这里只计算了中间Width范围内的数据,当不是第一行时,如下图左边所示,新的累加值只需减去移出的哪一行像素值同时加上移入的一行像素值,详见90到96
行代码。

  上面只计算了中间Width范围内的累加值,为了方便后续代码的编写以及使用SSE优化,下面的FillLeftAndRight_Mirror_C主要作用就是填充左边和右边分别填充数据,而且是按照镜像的方式进行填充。

  自从何凯明提出导向滤波后,因为其算法的不难性和管事,该算法获得了科普的应用,以至于新版的matlab都将其作为专业自带的函数之一了,利用他可以化解的享有的保边滤波器的能一举成功的难点,比如细节增强、HDR压缩、细节羽化、去雾、风格化,而且由于其保边特性,如果过多价值观函数中使用高斯滤波或者均值滤波的地点用他取代,能很好解决部分强边缘的连结不自然难题,比如retinex、Highlight/shadow等利用中,由此,神速的落到实处该算法具有很强的适用意义。

 

  还说一点,现在多数的CPU都扶助AVX256了,仍是可以选取AVX进一步加速,就好像代码该起来也不是很难,有趣味的意中人可以友善尝试。

 for (int Y = 0; Y < Height; Y++)
 {
     float PosY = Y * (SmallH - 1.0f) / Height;
     int YY = (int)PosY;
     float PartYY = PosY - YY;
     float InvertYY = 1 - PartYY;
     byte *LinePS = Guide + Y * Stride;
     byte *LinePD = Dest + Y * Stride;
     float *PtTopA = TempA + YY * Width * Channel;
     float *PtBottomA = PtTopA + Width * Channel;
     float *PtTopB = TempB + YY * Width * Channel;
     float *PtBottomB = PtTopB + Width * Channel;
     for (int X = 0;; X < Width * Channel; X++)
     {
         float ValueA = PtTopA[X] * InvertYY + PtBottomA[X] * PartYY;
         float ValueB = PtTopB[X] * InvertYY + PtBottomB[X] * PartYY;
         LinePD[X] = IM_ClampFHtoByte(ValueA * LinePS[X] + ValueB * 255);
     }
 }

     

  FillLeftAndRight_Mirror_C紧若是用来赢得两边镜像数据的(直接得到,不调用IM_GetMirrorPos函数),比如比如Array原始数据为
***abcdefgh*** (Length = 8, Radius =
3),则结果为dcbabcdefghgfe。

     
那样改动后,所有的boxfilter均是对下取样后的数码开展处理,当s=4时,计算量减弱到原始的1/16,而s=5,则裁减到了原始的1/25,当然那么些时候多了2个下取样和2个上取样的算法,下取样由于是压缩,统计量很小,无需关切,而上采样,计算量和原图大小有关,按照我的测评,那么些上采样的耗时可能占全体进程的形似时间左右的,是非凡值得注意优化的。

 for (int Y = 0; Y < Height; Y++)
 {
     float PosY = Y * (SmallH - 1.0f) / Height;
     int YY = (int)PosY;
     float PartYY = PosY - YY;
     float InvertYY = 1 - PartYY;
     byte *LinePS = Guide + Y * Stride;
     byte *LinePD = Dest + Y * Stride;
     float *PtTopA = TempA + YY * Width * Channel;
     float *PtBottomA = PtTopA + Width * Channel;
     float *PtTopB = TempB + YY * Width * Channel;
     float *PtBottomB = PtTopB + Width * Channel;
     for (int X = 0;; X < Width * Channel; X++)
     {
         float ValueA = PtTopA[X] * InvertYY + PtBottomA[X] * PartYY;
         float ValueB = PtTopB[X] * InvertYY + PtBottomB[X] * PartYY;
         LinePD[X] = IM_ClampFHtoByte(ValueA * LinePS[X] + ValueB * 255);
     }
 }

 

     

ca88亚洲城网站 9

    A3+P3+P2+P1+P0  A3+P2+P1+P0  A3+P1+P0  A3+P0;

ca88亚洲城网站 10

       6: q = meana. * +
meanb

  使用多个__m128i
变量首倘若为了丰盛利用XMM寄存器,其中SumI32函数如下,首假若为着总括__m128i
内七个整数的累加值。

       6: q = meana. * +
meanb

   
 关于上述浮点版本的Boxfilter,其实还有一种更好的落到实处情势。我在13行代码完结最便捷最便捷的积分图像算法中也提供了一段落成方框模糊的代码,当然分外代码还不是最优的,因为其中的pixlecount需求各种像素都重新总计,其实当半径较小时中间部分的像素的pixlecount为固定值,由此得以把边缘部分的像素特殊处理,对于本例,是急需举行的浮点版本的算法,那对于中等有些的
/ pixlecount操作就应该可以改为 *Invpixlecount,其中Invpixlecount =
1.0f/pixlecount,变除法为乘法,而且那部分盘算还是可以够很容易的用SSE落成。我测试过,立异后的贯彻和解析opencv中Box
Filter的完毕并指出越来越加速的方案(源码共享)
 
那篇文好章的速度方驾齐驱,但此间有个优势就是足以并行的。此外,最重点的少数时,当要总计上述浮点版半径版本的boxfilter时,积分图是不须求再一次重新统计的,而积分图的计算所占的耗时起码有一半左右。由此,那么些场面使用积分图版本的盒子滤波会更有优势。

  

ca88亚洲城网站 11

  注意最终的IM_ClampFHtoByte函数是将括号内的值限制在0和255以内的。

  镜像就是倒序的进度,直接运用SSE的shuffle函数很方便完成。

     

      在何的舆论中已经证实下采样比例 s
取4时,总结的结果和高精度结果也依然这么些靠近的,我在本人的完成里s
取到了5。

  总括数组的丰裕和优化也利于。

 
   我刚刚提的在去雾中自己实用的小Trick实际上就是第六步及第七步分歧,我的点子可发挥如下:

ca88亚洲城网站 12

ca88亚洲城网站 13

 for (int Y = 0; Y < SmallH; Y++)
 {
     float PosX = 0;
     float AddX = (SmallW - 1.0f) / Width;        //  主要是为了减少下面插值向右增1的指针超过范围,但这样做其实是和精确的算法有一点点差异的
     float *LinePDA = TempA + Y * Width * Channel;    //  TempA和TempB为临时分配的大小为(SmallH * Width * Channel * sizeof(float)大小的内存
     float *LinePDB = TempB + Y * Width * Channel;
     float *LinePA = MeanA + Y * SmallW * Channel;
     float *LinePB = MeanB + Y * SmallW * Channel;
     if (Channel == 1)
     {
         for (int X = 0; X < Width; X++)
         {
             int XX = (int)PosX;
             float PartXX = PosX - XX;
             float InvertXX = 1 - PartXX;
             float *PtLeftA = LinePA + XX;
             float *PtLeftB = LinePB + XX;
             LinePDA[X] = PtLeftA[0] * InvertXX + PtLeftA[1] * PartXX;
             LinePDB[X] = PtLeftB[0] * InvertXX + PtLeftB[1] * PartXX;
             PosX += AddX;
         }
     }
   //  ...................
 }

      何凯明在2015又公布了一篇《法斯特Guided Filter》的篇章,解说了一种很实用的更便捷的导向滤波流程:

 

     
很明显,由于I的插足统计,何的做法能更大程度上维持结果和原汁原味的近乎,而我的点子则会发出较大的块状相似,所以住户大神就是大神。

 

  彩色图工程:https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

   

       
 ca88亚洲城网站 14 
 ———–——>   ca88亚洲城网站 15 
 ———–——>  ca88亚洲城网站 16

  简单的叙述下各函数的功用先。

     
何凯明在导向滤波一文的相干资料中提供了其matlab代码,或者用上面的流程也可以清楚的发布:

     
有不少情人或者不亮堂,如果把地点的IM_ClampFHtoByte那个函数去掉,间接选拔括号内的代码,VS的编译器可以很好的对地方代码进行向量化编译(VS编译只要你未曾把代码生成–》启用增强指令集设置成无增强指令/arch:IA32,哪怕设置为未安装,都会把浮点的代码编译为SIMD相关指令的),而一旦大家对分化的Channel,比如3通道4坦途在循环里举行后,很糟糕,根据大家的进展循环的辩论,速度应该加紧,但实际却反倒了。所以大家要求丰裕明白编译器的向量化特性,就能写成更迅捷的代码。

int SumofArray_SSE(int *Array, int Length)
{
    int BlockSize = 8, Block = Length / BlockSize;
    __m128i Sum1 = _mm_setzero_si128();
    __m128i Sum2 = _mm_setzero_si128();
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        Sum1 = _mm_add_epi32(Sum1, _mm_loadu_si128((__m128i *)(Array + X + 0)));
        Sum2 = _mm_add_epi32(Sum2, _mm_loadu_si128((__m128i *)(Array + X + 4)));
    }
    int Sum = SumI32(_mm_add_epi32(Sum1, Sum2));
    //  处理剩余不能被SSE优化的数据
    return Sum;
}

  大家来看了地点的6次取mean计算的进程,也就是浮点数的boxfilter,那些事物已经是老掉牙的一个算法了,我在几年前探究过opencv内部的那个算法,并且提议了一种比opencv落成更快的法子,详见解析opencv中Box
Filter的完结并提议越来越加速的方案(源码共享)
 一文。不过那里的拍卖时针对字节数据的,其内部用到了一些整形数据的SSE优化,假设原本数据是浮点数,那反而就更为简单了,因为SSE指令生来就是为浮点数服务的。

ca88亚洲城网站 17

  最终就是100到104行的代码的变换。

     
使用SSE优化能将上述进程提速2倍以上。

   
 普通意义的boxfilter肯定是心有余而力不足支撑浮点半径的(那分化于高斯模糊),一种变化的主意就是取浮点半径前后的多个整形半径值做模糊,然后再线性插值,举个例子,假诺下取样后的半径为4.4,则分级总结R1
= boxfilter(4)以及R2 = boxfilter(5),最终合成得到结果R:

       
在更新了上述累加值后,大家早先拍卖总括均值了,对于每行的率先个点,SumofArray_C总括了前2*R
+
1个列的累加值,这么些累加值就是该点周边(2*R+1)*(2*R+1)个因素的累积和,对于一行的其它像素,其实就如同于行方向列累加值的立异,减去移出的加盟新进的列,如下图右边图所示,102到104行即落到实处了该进程。

   
 共享下一个C#做的Demo,以便供有趣味的情西洋参考比较: http://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

  那段代码用SSE去优化的祸害的脑细胞有点多,而且由于其统计量不是太大,意义可能有数。

  源代码下载:https://files.cnblogs.com/files/Imageshop/FastBlur.rar

  那段代码用SSE去优化的侵害的脑细胞有点多,而且由于其总括量不是太大,意义或许有限。

   
 共享下一个C#做的Demo,以便供有趣味的心上高丽参考相比较: http://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

void FillLeftAndRight_Mirror_SSE(int *Array, int Length, int Radius)
{
    int BlockSize = 4, Block = Radius / BlockSize;
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        __m128i SrcV1 = _mm_loadu_si128((__m128i *)(Array + Radius + Radius - X - 3));
        __m128i SrcV2 = _mm_loadu_si128((__m128i *)(Array + Radius + Length - X - 5));
        _mm_storeu_si128((__m128i *)(Array + X), _mm_shuffle_epi32(SrcV1, _MM_SHUFFLE(0, 1, 2, 3)));
        _mm_storeu_si128((__m128i *)(Array + Radius + Length + X), _mm_shuffle_epi32(SrcV2, _MM_SHUFFLE(0, 1, 2, 3)));
    }
    //  处理剩余不能被SSE优化的数据
}

     
不过即使是这样,由于6次计算以及中等的别样一些浮点运算,照旧给所有算法带来了很大的运算开支和内存开支,在无数场馆仍旧无法满足需要的,比如实时去雾等现象。在中期我的飞速去雾完成中,都是先拔取下采样图的导向滤波结果,然后再双线性插值放大得到大图的透射率图,即便在视觉效果上能一蹴即至去雾算法的进程难题,不过一旦是任何场景的导向滤波须求,仍旧会合到不少瑕疵的。

ca88亚洲城网站 18

int BlockSize = 4, Block = (Width - 1) / BlockSize;
__m128i OldSum = _mm_set1_epi32(LastSum);
__m128 Inv128 = _mm_set1_ps(Inv);
for (int X = 1; X < Block * BlockSize + 1; X += BlockSize)
{
    __m128i ColValueOut = _mm_loadu_si128((__m128i *)(ColValue + X - 1));
    __m128i ColValueIn = _mm_loadu_si128((__m128i *)(ColValue + X + Radius + Radius));
    __m128i ColValueDiff = _mm_sub_epi32(ColValueIn, ColValueOut);                            //    P3 P2 P1 P0                                                
    __m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4));        //    P3+P2 P2+P1 P1+P0 P0
    __m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8));                 //    P3+P2+P1+P0 P2+P1+P0 P1+P0 P0
    __m128i NewSum = _mm_add_epi32(OldSum, Value);
    OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3));                              //    重新赋值为最新值
    __m128 Mean = _mm_mul_ps(_mm_cvtepi32_ps(NewSum), Inv128);
    _mm_storesi128_4char(_mm_cvtps_epi32(Mean), LinePD + X);
}

 

 
   我正好提的在去雾中自己实用的小Trick实际上就是第六步及第七步不相同,我的章程可发挥如下:

  IM_GetMirrorPos函数紧倘若收获某一个岗位Pos按照镜像的格局处理时在Length动向的坐标,那里Pos可以为负值,那个至关重借使用来收获前期的坐标偏移的。      

     
 由于那些上采样是针对性浮点型的多少,所以中级的精度损失难题可以不要考虑,而假设是图像的字节数据,则要慎重了。

     
 由地点的率先个图到第一个图的光景代码如下:

 

__m128i SrcI = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i const *)(LinePS + X)), Zero);        //    Load the lower 64 bits of the value pointed to by p into the lower 64 bits of the result, zeroing the upper 64 bits of the result.
__m128 SrcFL = _mm_cvtepi32_ps(_mm_unpacklo_epi16(SrcI, Zero));                                //    转换为浮点
__m128 SrcFH = _mm_cvtepi32_ps(_mm_unpackhi_epi16(SrcI, Zero));

  我们见到了下边的6次取mean计算的进度,也就是浮点数的boxfilter,这几个东西已经是老掉牙的一个算法了,我在几年前切磋过opencv内部的那个算法,并且提议了一种比opencv完成更快的方法,详见ca88亚洲城网站,解析opencv中Box
Filter的完成并指出越来越加快的方案(源码共享)
 一文。可是那里的处理时针对字节数据的,其内部用到了一部分整形数据的SSE优化,固然原本数据是浮点数,那反而就尤其简单了,因为SSE指令生来就是为浮点数服务的。

  这一个时候3000*2000的图能做到25ms,,基本上接近自己革新的OPENCV的代码的快慢了。

  而由第三个图到第多个图的长河大概可有效下边的代码表述:

   
 在内存占用方面,也足以做大批量的优化工作,哪怕是对下取样图举办处理,第一,导向前必须把图像的字节数据归一化为0到1之间的浮点数据,很显眼,大家只要下采样大小的归一化数据,那么这一个进度就应当很自然的一向由原有大小图像投射到下取样的浮点数据,而毫无再在当中转来转去, 那么些下采样的内存占用大小为(W
* H )/(S * S) * channel * sizeof(float)
.第二,导向的高中级的各进程用到了汪洋的高中级变量,像原小编使用matlab的代码为了参数算法清楚,就是为每个中间数据分配内存,可是实际操作中,为节省资源,必须加以优化,我们注意阅览,就会发现有点变量用完就不会再也利用了,当导向图和原图差别时,我总括了只须要4
* (W * H )/(S * S) * channel *
sizeof(float)大小的内存,即使导向图和原图相同,则只须求2 * (W * H )/(S
* S) * channel *
sizeof(float),那些数目依旧含有下采样图的内存占用的吧。考虑在均值滤波里还必要一份附加大小的内存,以及最终混合时的为了涨价的
2 * (H / S) * W * channel *
sizeof(float)的内存,当S=4时加起来也就是原图多一点点的内存。

  SumofArray_C重如若测算一个数组的兼具的元素的和。

   
 如此处理后,在多数气象下(除了下取样后的半径为整数,比如原本半径为12,s=4,那是r’=3),统计量又会有些增添一些,须求统计小图的12次boxfilter了,不过何必纠结这些了吗。

   
 首先,第一,步骤6中的多个采样进度不要分开写,直接写到同一个for循环内部,那样可以省去成千成万坐标的测算进度,第二,那里一般的上采样平日选取双线性插值就OK了,网络上有很多有关双线性插值的SSE优化的代码,可是那些基本都是对准32位的图像做的优化,搬到24位和8位中是不适用的,而我们会在50%上述的票房价值中相见24位图像,所以说啊,网络上的事物虽多,但精华太少。

__m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4)); 这句的_mm_slli_si128得到中间结果 P2 P1 P0 0, 再和P3 P2 P1 P0相加得到

    Value_Temp =  P3+P2   P2+P1  P1+P0   P0

__m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8));这句的_mm_slli_si128得到中间结果P1+P0 P0 0 0,再和P3+P2 P2+P1 P1+P0  P0相加得到:

    Value = P3+P2+P1+P0   P2+P1+P0   P1+P0   P0
好简单的过程。

  OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3)); 这一句为什么要这样写,恐怕还是读者自己体会比较好,似乎不好用文字表达。

   
 关于上述浮点版本的Boxfilter,其实还有一种更好的贯彻方式。我在13行代码完成最连忙最急速的积分图像算法中也提供了一段已毕方框模糊的代码,当然卓殊代码还不是最优的,因为其中的pixlecount须要各种像素都重新计算,其实当半径较小时中间有些的像素的pixlecount为固定值,因而可以把边缘部分的像素特殊处理,对于本例,是内需展开的浮点版本的算法,那对于中等有些的
/ pixlecount操作就应当可以成为 *Invpixlecount,其中Invpixlecount =
1.0f/pixlecount,变除法为乘法,而且那有的盘算还足以很简单的用SSE完毕。我测试过,革新后的完毕和解析opencv中Box
Filter的贯彻并建议进一步加速的方案(源码共享)
 
那篇文好章的速度齐镳并驱,但那里有个优势就是足以相互的。别的,最关键的一些时,当要总结上述浮点版半径版本的boxfilter时,积分图是不须要再一次重复计算的,而积分图的统计所占的耗时起码有一半左右。因而,那几个场面使用积分图版本的盒子滤波会更有优势。

__m128i SrcI = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i const *)(LinePS + X)), Zero);        //    Load the lower 64 bits of the value pointed to by p into the lower 64 bits of the result, zeroing the upper 64 bits of the result.
__m128 SrcFL = _mm_cvtepi32_ps(_mm_unpacklo_epi16(SrcI, Zero));                                //    转换为浮点
__m128 SrcFH = _mm_cvtepi32_ps(_mm_unpackhi_epi16(SrcI, Zero));

  __m128i ColValueDiff =
_mm_sub_epi32(ColValueIn, ColValueOut);
这句代码得到了移入和移出种类的差值,大家计为;  P3 P2 P1 P0
(高位到低位)由于算法的增加特性,若是说OldSum的原始值为A3 A3 A3
A3,那么新的NewSum就应当是:

   
 在内存占用方面,也足以做大批量的优化办事,哪怕是对下取样图进行拍卖,第一,导向前必须把图像的字节数据归一化为0到1里面的浮点数据,很明显,大家即使下采样大小的归一化数据,那么这些进度就活该很当然的直接由原来大小图像投射到下取样的浮点数据,而不用再在中间转来转去, 这几个下采样的内存占用大小为(W
* H )/(S * S) * channel * sizeof(float)
.第二,导向的中档的各进程用到了汪洋的中级变量,像原作者使用matlab的代码为了参数算法清楚,就是为每个中间数据分配内存,不过实际操作中,为节省资源,必须加以优化,大家注意寓目,就会发觉有点变量用完就不会再一次利用了,当导向图和原图不一致时,我总括了只要求4
* (W * H )/(S * S) * channel *
sizeof(float)大小的内存,若是导向图和原图相同,则只须求2 * (W * H )/(S
* S) * channel *
sizeof(float),这么些数据或者含有下采样图的内存占用的呢。考虑在均值滤波里还亟需一份附加大小的内存,以及尾声混合时的为了涨价的
2 * (H / S) * W * channel *
sizeof(float)的内存,当S=4时加起来也就是原图多一点点的内存。

   
 由于在盘算进度中真的存在一些结出出乎了0和255的界定,因而如果把IM_ClampFHtoByte函数去除,对有些图像会油可是生噪点,由此,大家不可以完全依靠编译器的向量化优化了,那就非得自己写SIMD指令,由于SIMD自带了饱和处理的连带函数,而上述内部的X
的for循环是很不难用SSE处理的,唯一必要小心的就是要求把LinePS对应的字节数据转换为浮点数据,那里自己概括的提示可以用如下指令将8个字节数据转换为8个浮点数:

  能够说,近日那一个速度已经基本上达到了CPU的巅峰了,不过测试过IPP的速度,就如比这一个还要快点,不消除他利用了AVX,也不免除他选择多核的资源。

       7:   q = fupsample(q,
s)

   
 如此处理后,在大部分处境下(除了下取样后的半径为整数,比如原本半径为12,s=4,那是r’=3),总结量又会稍为增添某些,须要总计小图的12次boxfilter了,可是何必纠结这几个了吧。

  不佳意思,图太小,速度为0ms了。

      在何的诗歌中曾经证实下采样比例 s
取4时,总计的结果和纯粹结果也仍然十分靠近的,我在自家的落实里s
取到了5。

     
那里如此做的此外一个便宜是在Y循环中计算是独自的,由此都可以使用OPENMP加快。

  1 /// <summary>
  2 /// 按照Tile方式进行数据的扩展,得到扩展后在原尺寸中的位置,以0为下标。
  3 /// </summary>
  4 int IM_GetMirrorPos(int Length, int Pos)
  5 {
  6     if (Pos < 0)
  7         return -Pos;
  8     else if (Pos >= Length)
  9         return Length + Length - Pos - 2;
 10     else    
 11         return Pos;
 12 }
 13 
 14 void FillLeftAndRight_Mirror_C(int *Array, int Length, int Radius)
 15 {
 16     for (int X = 0; X < Radius; X++)
 17     {
 18         Array[X] = Array[Radius + Radius - X];
 19         Array[Radius + Length + X] = Array[Radius + Length - X - 2];
 20     }
 21 }
 22 
 23 int SumofArray_C(int *Array, int Length)
 24 {
 25     int Sum = 0;
 26     for (int X = 0; X < Length; X++)
 27     {
 28         Sum += Array[X];
 29     }
 30     return Sum;
 31 }
 32 
 33 /// <summary>
 34 /// 将整数限幅到字节数据类型。
 35 /// </summary>
 36 inline unsigned char IM_ClampToByte(int Value)            //    现代PC还是这样直接写快些
 37 {
 38     if (Value < 0)
 39         return 0;
 40     else if (Value > 255)
 41         return 255;
 42     else
 43         return (unsigned char)Value;
 44     //return ((Value | ((signed int)(255 - Value) >> 31)) & ~((signed int)Value >> 31));
 45 }
 46 
 47 int IM_BoxBlur_C(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
 48 {
 49     int Channel = Stride / Width;
 50     if ((Src == NULL) || (Dest == NULL))                    return IM_STATUS_NULLREFRENCE;
 51     if ((Width <= 0) || (Height <= 0) || (Radius <= 0))        return IM_STATUS_INVALIDPARAMETER;
 52     if (Radius < 1)                                            return IM_STATUS_INVALIDPARAMETER;
 53     if ((Channel != 1) && (Channel != 3) && (Channel != 4))    return IM_STATUS_NOTSUPPORTED;
 54 
 55     Radius = IM_Min(IM_Min(Radius, Width - 1), Height - 1);        //    由于镜像的需求,要求半径不能大于宽度或高度-1的数据
 56 
 57     int SampleAmount = (2 * Radius + 1) * (2 * Radius + 1);
 58     float Inv = 1.0 / SampleAmount;
 59 
 60     int *ColValue = (int *)malloc((Width + Radius + Radius) * (Channel == 1 ? Channel : 4) * sizeof(int));
 61     int *ColOffset = (int *)malloc((Height + Radius + Radius) * sizeof(int));
 62     if ((ColValue == NULL) || (ColOffset == NULL))
 63     {
 64         if (ColValue != NULL)    free(ColValue);
 65         if (ColOffset != NULL)    free(ColOffset);
 66         return IM_STATUS_OUTOFMEMORY;
 67     }
 68     for (int Y = 0; Y < Height + Radius + Radius; Y++)
 69         ColOffset[Y] = IM_GetMirrorPos(Height, Y - Radius);
 70 
 71     if (Channel == 1)
 72     {
 73         for (int Y = 0; Y < Height; Y++)
 74         {
 75             unsigned char *LinePD = Dest + Y * Stride;
 76             if (Y == 0)
 77             {
 78                 memset(ColValue + Radius, 0, Width * sizeof(int));
 79                 for (int Z = -Radius; Z <= Radius; Z++)
 80                 {
 81                     unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride;
 82                     for (int X = 0; X < Width; X++)
 83                     {
 84                         ColValue[X + Radius] += LinePS[X];                                            //    更新列数据
 85                     }
 86                 }
 87             }
 88             else
 89             {
 90                 unsigned char *RowMoveOut = Src + ColOffset[Y - 1] * Stride;                //    即将减去的那一行的首地址    
 91                 unsigned char *RowMoveIn = Src + ColOffset[Y + Radius + Radius] * Stride;    //    即将加上的那一行的首地址
 92                 for (int X = 0; X < Width; X++)
 93                 {
 94                     ColValue[X + Radius] -= RowMoveOut[X] - RowMoveIn[X];                                            //    更新列数据
 95                 }
 96             }
 97             FillLeftAndRight_Mirror_C(ColValue, Width, Radius);                //    镜像填充左右数据
 98             int LastSum = SumofArray_C(ColValue, Radius * 2 + 1);                //    处理每行第一个数据                                
 99             LinePD[0] = IM_ClampToByte(LastSum * Inv);
100             for (int X = 1; X < Width; X++)
101             {
102                 int NewSum = LastSum - ColValue[X - 1] + ColValue[X + Radius + Radius];
103                 LinePD[X] = IM_ClampToByte(NewSum * Inv);
104                 LastSum = NewSum;
105             }
106         }
107     }
108     else if (Channel == 3)
109     {
110 
111     }
112     else if (Channel == 4)
113     {
114 
115     }
116     free(ColValue);
117     free(ColOffset);
118     return IM_STATUS_OK;
119 }

   
 普通意义的boxfilter肯定是无法支撑浮点半径的(那不一致于高斯模糊),一种变更的方法就是取浮点半径前后的七个整形半径值做模糊,然后再线性插值,举个例证,要是下取样后的半径为4.4,则分别计算R1
= boxfilter(4)以及R2 = boxfilter(5),最后合成获得结果R:

     
何凯明在导向滤波一文的连锁材料中提供了其matlab代码,或者用下边的流水线也得以清晰的抒发:

  以前没察觉到,就像此不难的代码用C写后能暴发速度也是很诱人的,3000*2000的图能不辱义务39ms,如若在编译选项里勾选编译器的“启用增强指令集:流式处理
SIMD 扩充 2 (/arch:SSE2)”,
则系统会对上述所有浮点总计的片段使用有关的SIMD指令优化,如下图所示:

  本文简要的笔录了自我在优化导向滤波完成的历程中所适用的优化措施和一些细节,以防时间久了后自己都不记得了,不过请不要向自家一直索取源代码。

               R = R1 * (1 – 0.4) + R2
* 0.4;

   最终一个_mm_storesi128_4char是自家要好定义的一个将1个__m128i里面的4个32位整数转换为字节数并保留到内存中的函数,详见附件代码。

ca88亚洲城网站 19

   
 在一台I5的记录本上,选取默许设置,以自家为导向图处理3000*2000的24位图像必要约55ms,假若是灰度图大概是20ms,那几个和优化后的
boxblur速度基本一致,借使开启动二十四线程,比如开七个线程,仍能提速25%左右,再多也无辅助了。

int BlockSize = 8, Block = Width / BlockSize;
__m128i Zero = _mm_setzero_si128();
for (int X = 0; X < Block * BlockSize; X += BlockSize)
{
    int *DestP = ColValue + X + Radius;
    __m128i MoveOut = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(RowMoveOut + X)), Zero);
    __m128i MoveIn = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(RowMoveIn + X)), Zero);
    __m128i Diff = _mm_sub_epi16(MoveIn, MoveOut);                        //    注意这个有负数也有正数的,有负数时转换为32位是不能用_mm_unpackxx_epi16体系的函数
    _mm_storeu_si128((__m128i *)DestP, _mm_add_epi32(_mm_loadu_si128((__m128i *)DestP), _mm_cvtepi16_epi32(Diff)));
    _mm_storeu_si128((__m128i *)(DestP + 4), _mm_add_epi32(_mm_loadu_si128((__m128i *)(DestP + 4)), _mm_cvtepi16_epi32(_mm_srli_si128(Diff, 8))));
}
//  处理剩余不能被SSE优化的数据

     
里面的浮点统计的进度的SSE代码就和平凡的函数调用没什么却别,最终的写到LinePD那一个字节数据的进度可以用_mm_storel_epi64以及关于活动搞定。

   

  代码不解释。

     
我使用的一个优化措施时,先举行水平方向的上采样到一个缓冲区中(Width  *
SmallH),然后在从那一个缓冲区中沿着中度方向缓冲到(Width *
Height),如下图所示:

     
别的一个标题,在上头的流程2的率先步中,对boxfilter的半径r也是开展了同比例的收缩的,注意到boxfilter的半径平常状态下我们都是用的平头,倘使缩短后的r’也进行取整的话,举例来说,对于s
=4的图景下,半径为8、9、10、11那多样境况最终获得的导向滤波结果就全盘等同了,如同那不符合大家对算法严酷性的必要,所以大家要援救一种浮点半径的boxfilter。

  IM_BoxBlur_C函数内部即为模糊的函数体,选用的优化思路完全和自由半径中值滤波(扩张至百分比滤波器)O(1)时间复杂度算法的规律、完毕及成效是同等的。当半径为R时,方框模糊的值等于以某点为主导,左右左右各伸张R像素的的范围内拥有像素的概括,像素总个数为(2*R+1)*(2*R+1)个,当然咱们也可以把她分成(2*R+1)列,每列有(2*R+1)个元素本例的优化措施我们就是把累加数据分为一列一列的,丰盛利用重复音讯来达成速度升高。

 

     
 由于那几个上采样是针对浮点型的数目,所以中间的精度损失难题可以毫无考虑,而一旦是图像的字节数据,则要慎重了。

  

//    Convert 4 32-bit integer values to 4 unsigned char values.
void _mm_storesi128_4char(__m128i Src, unsigned char *Dest)
{
    __m128i T = _mm_packs_epi32(Src, Src);
    T = _mm_packus_epi16(T, T);
    *((int*)Dest) = _mm_cvtsi128_si32(T);
}

相关文章