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

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

      在参考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结构的完整代码,我以后的文章都将以改结构为依托进行处理。

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

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

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

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

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏CDA数据分析师

工具 | R语言数据可视化之数据分布图(直方图、密度曲线、箱线图、等高线、2D密度图)

数据分布图简介 绘制基本直方图 基于分组的直方图 绘制密度曲线 绘制基本箱线图 往箱线图添加槽口和均值 绘制2D等高线 绘制2D密度图 数据分布图简介 中医上讲...

21910
来自专栏ml

nyoj-----幸运三角形

幸运三角形 时间限制:1000 ms  |  内存限制:65535 KB 难度:3 描述         话说有这么一个图形,只有两种符号组成(‘+’或者‘-’...

28310
来自专栏PPV课数据科学社区

【工具】SAS 常用函数汇总

? 一、数学函数 ABS(x) 求x的绝对值。 MAX(x1,x2,…,xn) 求所有自变量中的最大一个。 MIN(x1,x2,…,xn) 求所有自变量...

2273
来自专栏大数据挖掘DT机器学习

tensorflow LSTM + CTC实现端到端OCR

本文github源码地址: 在公众号 datadw 里 回复 OCR 即可获取。 最近在做OCR相关的东西,关于OCR真的是有悠久了历史了,最开始用tes...

6073
来自专栏大数据风控

R中时间序列分析-趋势分析Trend

趋势分析(Trend) 常用趋势的数学函数 线性函数 y=ax+b 指数函数 y=a^x 二次函数 y=ax^2+bx+c 曲线拟合方法 nls...

1838
来自专栏数说工作室

SAS里的平行世界 | 【SAS Says · 扩展篇】IML:1.入门

本节目录: 1. IML基本语句 2. 标量、向量与矩阵 (1)定义标量 (2)定义向量 (3)定义矩阵 3. 矩阵的元素运算 (1...

3476
来自专栏用户2442861的专栏

openCV—访问与操作像素(图片的区域的更改)

结果如下: 位置(0,0)处的像素 - 红:225,绿:138,蓝:128 位置(0,0)处的像素 - 红:200,绿:150,蓝:100

591
来自专栏闪电gogogo的专栏

IEEE Trans 2009 Stagewise Weak Gradient Pursuits论文学习

论文在第二部分先提出了贪婪算法框架,如下截图所示: ? 接着根据原子选择的方法不同,提出了SWOMP(分段弱正交匹配追踪)算法,以下部分为转载《压缩感知重构算法...

3208
来自专栏程序员叨叨叨

2.4 图形硬件

这一节中主要阐述图形硬件的相关知识,主要包括 GPU 中数据的存放硬件, 以及各类缓冲区的具体含义和用途,如:z buffer(深度缓冲区)、stencil b...

592
来自专栏落影的专栏

Metal图像处理——直方图均衡化

首先,我们用直方图来表示一张图像:横坐标代表的是颜色值,纵坐标代表的是该颜色值在图像中出现次数。

1105

扫码关注云+社区