文档章节

任意半径局部直方图类算法在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
粉丝 64
博文 6196
码字总数 1876
作品 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
OpenCV之LBP算法学习

目录 前言 LBP算法概述 指局部二值模式,英文全称:,是一种用来描述图像局部特征的算子,特征具有灰度不变性和旋转不变性等显著优点。常应用于人脸识别和目标检测中,在中有使用特征进行人脸...

simonforfuture
2017/12/20
0
0
LBP特征原理及代码实现

一、LBP特征的背景介绍 LBP指局部二值模式,英文全称:Local Binary Pattern,是一种用来描述图像局部特征的算子,LBP特征具有灰度不变性和旋转不变性等显著优点。它是由T. Ojala, M.Pietikä...

Quincuntial
2016/01/19
0
0
机器学习之聚类

聚类 1、聚类试图将数据集中的数据划分为若干个通常是不相交的子集,每个子集称为一个“簇”(cluster) 2、聚类过程仅能自动形成簇结构,簇所对应的概念语义需由使用者来把握和命名 思考:簇所...

qq_37634812
2017/12/01
0
0

没有更多内容

加载失败,请刷新页面

加载更多

jdk 1.8 在线文档

https://docs.oracle.com/javase/8/docs/api/

kuchawyz
7分钟前
1
0
python:有个叫strip的东西.....

strip.....脱.......呃,这个嘛,好吧,也许python的开发团队并不忌讳strip这个词的意思[] lstrip() >>> b' spacious '.lstrip()b'spacious '>>> b'www.example.com'.lstrip(b'cmowz......

Oh_really
20分钟前
0
0
Rails 用现代 Rails 逃离单页面应用 “兔子洞”

在工作共总是觉得turbolinks非常爽,但是却总是被说成是过时的技术,大家都喜欢spa,哪怕不用的spa的人也是禁用掉的多,找不到很好的理由劝说别人使用,这篇文章说的很到位,或者说至少是牛人...

wmzsonic
40分钟前
0
0
Hive 分布式搭建,Spark集成Hive记录

本帖详细介绍搭建步骤,仅仅记录自己搭建过程以及采坑经历。 前提环境: Hadoop集群 版本2.7.2 Spark集群 版本2.1.0 Linux版本 Centos7 准备搭建 MySql版本5.5.61 ,Hive-2.1.0 去官网下载M...

我爱春天的毛毛雨
43分钟前
3
0
打包QML程序

1、windeployqt执行路径(D:\Qt\5.12.0\msvc2017_64\bin)加入到PATH中 2、使用Qt自带的命令行交互 Command 终端(Qt 5.12.0 64-bit for Desktop (MSVC 2017))切换到 Release 编译成功的exe...

渣渣曦
今天
4
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部