文档章节

任意半径局部直方图类算法在PC中快速实现的框架。

abcijkxyz
 abcijkxyz
发布于 2016/11/22 16:39
字数 2018
阅读 1
收藏 0

 

     代码下载地址:http://files.cnblogs.com/files/Imageshop/BaseFile.rar

     在图像处理中,局部算法一般来说,在很大程度上会获得比全局算法更为好的效果,因为他考虑到了图像领域像素的信息,而很多局部算法可以借助于直方图获得加速。同时,一些常规的算法,比如中值滤波、最大值滤波、最小值滤波、表面模糊等等都可以通过局部直方图进行加速。而传统的获取局部直方图计算量很大,特别是半径增加时,耗时会成平方关系增加。一些局部算法只有在半径较大时才会获得很好的效果,因此,必须找到一种合适的加速计算局部直方图的方式。

      在参考Median Filter in Constant Time.pdf一文附带的C的代码的基础上,本文提出了基于SSE加速的恒长任意半径局部直方图获取技术,可以大大加速算法的计算时间,特别是大半径时的提速更为明显。

      主要的优化思路是,沿着列方向一行一行的更行整行的列直方图,新的一行对应的列直方图更新时只需要减去已经不再范围内的那个像素同时加入新进入的像素的直方图信息。之后,对于一行中的第一个像素点,累加半径辐射范围内的列直方图,得到改点的局部直方图,对于行中的其他的像素,则类似于更新行直方图,先减去不在范围内那列的列直方图,然后加上移入范围内的列直方图。由于采用了基于SSE函数的加速过程,直方图想加和相减的速度较普通的加减法有了10倍以上的提速,因此大大的提高了整体的实用性。

       具体的过程我用代码加以说明:

    1、一些公用的内存分配过程

    TMatrix *Row = NULL, *Col = NULL; unsigned char *LinePS, *LinePD; int  X, Y, K, Width = Src->Width, Height = Src->Height; int *RowOffset, *ColOffSet; unsigned short *ColHist    = (unsigned short *)IS_AllocMemory(256 * (Width + 2 * Radius) * sizeof(unsigned short), true); if (ColHist == NULL) {Ret = IS_RET_ERR_OUTOFMEMORY; goto Done8;} unsigned short *Hist    = (unsigned short *)IS_AllocMemory(256 * sizeof(unsigned short), true); if (Hist == NULL) {Ret = IS_RET_ERR_OUTOFMEMORY; goto Done8;} Ret = GetValidCoordinate(Width, Height, Radius, Radius, Radius, Radius, Edge, &Row, &Col);        // 获取坐标偏移量
    if (Ret != IS_RET_OK) goto Done8;

  其中的ColHist用于保存一行像素对应的列直方图 ,注意这里的行是用的扩展后的行的大小即:Width + 2 * Radius。IS_AllocMemory是个内部使用了_mm_malloc定义的内存分配函数,主要是考虑SSE函数的16字节对齐问题。

      Hist变量用于保存每个像素点的局部直方图数据,任何基于局部直方图技术的函数最终都演变为对于该函数进行各种各样的计算。    

      GetValidCoordinate是一个用于辅助边界处像素点处理的函数,具体可详见附件中给出的代码。

      2、更新一行像素的列直方图

for (Y = 0; Y < Height; Y++) { if (Y == 0)                                            // 第一行的列直方图,要重头计算
 { for (K = -Radius; K <= Radius; K++) { LinePS = Src->Data + ColOffSet[K] * Src->WidthStep; for (X = -Radius; X < Width + Radius; X++) { ColHist[X * 256 + LinePS[RowOffset[X]]]++; } } } else                                                // 其他行的列直方图,更新就可以了
 { LinePS = Src->Data + ColOffSet[Y - Radius - 1] * Src->WidthStep; for (X = -Radius; X < Width + Radius; X++)        // 删除移出范围内的那一行的直方图数据
 { ColHist[X * 256 + LinePS[RowOffset[X]]]--; } LinePS = Src->Data + ColOffSet[Y + Radius] * Src->WidthStep; for (X = -Radius; X < Width + Radius; X++)        // 增加进入范围内的那一行的直方图数据
 { ColHist[X * 256 + LinePS[RowOffset[X]]]++; } }
  //  依次获取一行每个像素的局部直方图
//  根据局部直方图获的结果
}

  可见,这部分和普通的局部优化方式类似,没有什么特殊的地方。

  3、依次获取一行每个像素的局部直方图

    for (Y = 0; Y < Height; Y++) {      //  更新一行像素的列直方图 memset(Hist, 0, 256 * sizeof(unsigned short));        // 每一行直方图数据清零先 LinePS = Src->Data + Y * Src->WidthStep; LinePD = Dest->Data + Y * Dest->WidthStep; for (X = 0; X < Width; X++) { if (X == 0) { for (K = -Radius; K <= Radius; K++)            // 行第一个像素,需要重新计算 
                    HistgramAddShort(ColHist + K * 256, Hist); } else {   /*  HistgramAddShort(ColHist + RowOffset[X + Radius] * 256, Hist);   HistgramSubShort(ColHist + RowOffset[X - Radius - 1] * 256, Hist);        */ HistgramSubAddShort(ColHist + RowOffset[X - Radius - 1] * 256, ColHist + RowOffset[X + Radius] * 256, Hist);  // 行内其他像素,依次删除和增加就可以了
 }
        //  根据局部直方图获的结果
 LinePS++; LinePD++; } }

  上面处理的过程其实和2的过程的优化道理是类似的,只不过一个是行方向,一个是列方向,聪明者自然能明白,稍微愚钝者请自己多多斟酌,自然有豁然开朗的时刻。

  4、 根据局部直方图获的结果

  根据不同的算法需求,结合局部直方图信息来获取结果,比如最大值算法可以用如下方式获得:

for (K = 255; K >= 0; K--)
    {
        if (Hist[K] != 0)
        {
            LinePD[X] = K;
            break;
        }
    }

     关于直方图累加的代码如下:

/// <summary>
/// 无符号短整形直方图数据相加,Y = X + Y, 整理时间2014.12.28; 
/// </summary>
/// <param name="X">加数。</param>
/// <param name="Y">被加数,结果保存于该数中。</param>
/// <remarks>使用了SSE优化。</remarks>
void HistgramAddShort(unsigned short *X, unsigned short *Y)
{
    *(__m128i*)(Y + 0)        =    _mm_add_epi16(*(__m128i*)&Y[0],        *(__m128i*)&X[0]);        //    不要想着用自己写的汇编超过他的速度了,已经试过了
    *(__m128i*)(Y + 8)        =    _mm_add_epi16(*(__m128i*)&Y[8],        *(__m128i*)&X[8]);
    *(__m128i*)(Y + 16)        =    _mm_add_epi16(*(__m128i*)&Y[16],    *(__m128i*)&X[16]);
    *(__m128i*)(Y + 24)        =    _mm_add_epi16(*(__m128i*)&Y[24],    *(__m128i*)&X[24]);
    *(__m128i*)(Y + 32)        =    _mm_add_epi16(*(__m128i*)&Y[32],    *(__m128i*)&X[32]);
    *(__m128i*)(Y + 40)        =    _mm_add_epi16(*(__m128i*)&Y[40],    *(__m128i*)&X[40]);
    *(__m128i*)(Y + 48)        =    _mm_add_epi16(*(__m128i*)&Y[48],    *(__m128i*)&X[48]);
    *(__m128i*)(Y + 56)        =    _mm_add_epi16(*(__m128i*)&Y[56],    *(__m128i*)&X[56]);
    *(__m128i*)(Y + 64)        =    _mm_add_epi16(*(__m128i*)&Y[64],    *(__m128i*)&X[64]);
    *(__m128i*)(Y + 72)        =    _mm_add_epi16(*(__m128i*)&Y[72],    *(__m128i*)&X[72]);
    *(__m128i*)(Y + 80)        =    _mm_add_epi16(*(__m128i*)&Y[80],    *(__m128i*)&X[80]);
    *(__m128i*)(Y + 88)        =    _mm_add_epi16(*(__m128i*)&Y[88],    *(__m128i*)&X[88]);
    *(__m128i*)(Y + 96)        =    _mm_add_epi16(*(__m128i*)&Y[96],    *(__m128i*)&X[96]);    
    *(__m128i*)(Y + 104)    =    _mm_add_epi16(*(__m128i*)&Y[104],    *(__m128i*)&X[104]);
    *(__m128i*)(Y + 112)    =    _mm_add_epi16(*(__m128i*)&Y[112],    *(__m128i*)&X[112]);
    *(__m128i*)(Y + 120)    =    _mm_add_epi16(*(__m128i*)&Y[120],    *(__m128i*)&X[120]);
    *(__m128i*)(Y + 128)    =    _mm_add_epi16(*(__m128i*)&Y[128],    *(__m128i*)&X[128]);
    *(__m128i*)(Y + 136)    =    _mm_add_epi16(*(__m128i*)&Y[136],    *(__m128i*)&X[136]);
    *(__m128i*)(Y + 144)    =    _mm_add_epi16(*(__m128i*)&Y[144],    *(__m128i*)&X[144]);
    *(__m128i*)(Y + 152)    =    _mm_add_epi16(*(__m128i*)&Y[152],    *(__m128i*)&X[152]);
    *(__m128i*)(Y + 160)    =    _mm_add_epi16(*(__m128i*)&Y[160],    *(__m128i*)&X[160]);
    *(__m128i*)(Y + 168)    =    _mm_add_epi16(*(__m128i*)&Y[168],    *(__m128i*)&X[168]);
    *(__m128i*)(Y + 176)    =    _mm_add_epi16(*(__m128i*)&Y[176],    *(__m128i*)&X[176]);
    *(__m128i*)(Y + 184)    =    _mm_add_epi16(*(__m128i*)&Y[184],    *(__m128i*)&X[184]);
    *(__m128i*)(Y + 192)    =    _mm_add_epi16(*(__m128i*)&Y[192],    *(__m128i*)&X[192]);
    *(__m128i*)(Y + 200)    =    _mm_add_epi16(*(__m128i*)&Y[200],    *(__m128i*)&X[200]);
    *(__m128i*)(Y + 208)    =    _mm_add_epi16(*(__m128i*)&Y[208],    *(__m128i*)&X[208]);
    *(__m128i*)(Y + 216)    =    _mm_add_epi16(*(__m128i*)&Y[216],    *(__m128i*)&X[216]);
    *(__m128i*)(Y + 224)    =    _mm_add_epi16(*(__m128i*)&Y[224],    *(__m128i*)&X[224]);    
    *(__m128i*)(Y + 232)    =    _mm_add_epi16(*(__m128i*)&Y[232],    *(__m128i*)&X[232]);
    *(__m128i*)(Y + 240)    =    _mm_add_epi16(*(__m128i*)&Y[240],    *(__m128i*)&X[240]);
    *(__m128i*)(Y + 248)    =    _mm_add_epi16(*(__m128i*)&Y[248],    *(__m128i*)&X[248]);
}

  _mm_add_epi16可以一次性完成16个short类型的数据的加法,比传统的add指令快了很多倍。

     由于_mm_add_epi16是这对短整形数据进行的处理,因此,一般情况下改指令所能处理的半径不能大于127,如果需要大于127,则需要修改过程序中的short类型为int,同时需要使用_mm_add_epi32指令,这样程序的速度会有所下降。

  经过测试,在我的I5的台式机中,1024*768图像在直方图更新上所需要的平均之间约为30ms,相比局部算法的核心就算部分时间(比如上述的求最大值),可能大部分耗时并不在这里。

     附件的代码中有个完整的测试工程,并有我目前所有的TMatrix结构的完整代码,我以后的文章都将以改结构为依托进行处理。

     代码还共享了很多处理的函数,我很自信一定值得朋友去学习的。

     这种前后依赖的算法都有一个很致命的缺点,就是不可以并行,把图像分段处理,也会造成过多初始化耗时。

     

 

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

 

本文转载自:http://www.cnblogs.com/Imageshop/p/4441972.html

共有 人打赏支持
abcijkxyz
粉丝 63
博文 6196
码字总数 1876
作品 0
深圳
项目经理
私信 提问
CVPR论文《100+ Times Faster Weighted Median Filter (WMF)》的实现和解析(附源代码)。

  四年前第一次看到《100+ Times FasterWeighted Median Filter (WMF)》一文时,因为他附带了源代码,而且还是CVPR论文,因此,当时也对代码进行了一定的整理和解读,但是当时觉得这个算法...

Imageshop
2018/11/13
0
0
【火炉炼AI】机器学习055-使用LBP直方图建立人脸识别器

【火炉炼AI】机器学习055-使用LBP直方图建立人脸识别器 (本文所使用的Python库和版本号: Python 3.6, Numpy 1.14, scikit-learn 0.19, matplotlib 2.2 ) 在我前面的博文【火炉炼AI】机器学习...

炼丹老顽童
2018/11/01
0
0
图像特征提取三大法宝:HOG特征,LBP特征,Haar特征

(一)HOG特征 1、HOG特征: 方向梯度直方图(Histogram of Oriented Gradient, HOG)特征是一种在计算机视觉和图像处理中用来进行物体检测的特征描述子。它通过计算和统计图像局部区域的梯度...

ranjiewen
2016/09/20
0
0
目标检测的图像特征提取之(二)LBP特征

原文地址:http://blog.csdn.net/zouxy09/article/details/7929531 LBP(Local Binary Pattern,局部二值模式)是一种用来描述图像局部纹理特征的算子;它具有旋转不变性和灰度不变性等显著的...

迈克老狼1
2013/11/22
0
0
机器学习实战——LBP特征提取

全篇概述: LBP(Local Binary Pattern)算法 是一种描述图像特征像素点与各个像素点之间的灰度关系的局部特征的非参数算法,同时也是一张高效的纹理描述算法。 纹理是物体表面的自然特性,它描...

机器学习算法工程师
2018/11/05
0
0

没有更多内容

加载失败,请刷新页面

加载更多

开始看《Java学习笔记》

虽然书买了很久,但一直没看。这其中也写过一些Java程序,但都是基于IDE的帮助和对C#的理解来写的,感觉不踏实。 林信良的书写得蛮好的,能够帮助打好基础,看得出作者是比较用心的。 第1章概...

max佩恩
昨天
11
0
Redux 三大原则

1.单一数据源 在传统的MVC架构中,我们可以根据需要创建无数个Model,而Model之间可以互相监听、触发事件甚至循环或嵌套触发事件,这些在Redux中都是不被允许的。 因为在Redux的思想里,一个...

wenxingjun
昨天
7
0
跟我学Spring Cloud(Finchley版)-12-微服务容错三板斧

至此,我们已实现服务发现、负载均衡,同时,使用Feign也实现了良好的远程调用——我们的代码是可读、可维护的。理论上,我们现在已经能构建一个不错的分布式应用了,但微服务之间是通过网络...

周立_ITMuch
昨天
4
0
XML

学习目标  能够说出XML的作用  能够编写XML文档声明  能够编写符合语法的XML  能够通过DTD约束编写XML文档  能够通过Schema约束编写XML文档  能够通过Dom4j解析XML文档 第1章 xm...

stars永恒
昨天
2
0
RabbitMQ学习(2)

1. 生产者客户端 void basicPublish(String exchange, String routingKey, boolean mandatory, boolean immediate, BasicProperties props, byte[] body) 1. 在生产者客户端发送消息时,首先......

江左煤郎
昨天
4
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部