SSE图像算法优化系列十二:多尺度的图像细节提升。

无意中浏览一篇文章,中间提到了基于多尺度的图像的细节提升算法,尝试了一下,还是有一定的效果的,结合最近一直研究的SSE优化,把算法的步骤和优化过程分享给大家。

论文的全名是DARK IMAGE ENHANCEMENT BASED ON PAIRWISE TARGET CONTRAST AND MULTI-SCALE DETAIL BOOSTING,好像在百度上搜索不到,由于博客的空间不多了,这里就不上传了, 我贴出论文核心的字段。

  论文的核心思想类似于Retinex,使用了三个尺度的高斯模糊,再和原图做减法,获得不同程度的细节信息,然后通过一定的组合方式把这些细节信息融合到原图中,从而得到加强原图信息的能力。

值得一提的就是对D1的系数做了特殊的处理,这个是值得学习的。

这个算法的编码实在是简单,一个简单的C语言代码如下:

int IM_MultiScaleSharpen(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                                return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                                    return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3) && (Channel != 4))                return IM_STATUS_INVALIDPARAMETER;

    int Status = IM_STATUS_OK;

    unsigned char *B1 = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));        //    如果最后的大图在这里内存有问题,一种解决办法就是下面的模糊用BoxBlur
    unsigned char *B2 = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));        //    然后不要调用现有的BoxBlur函数,而是直接在本函数实现三种不同半径的模糊
    unsigned char *B3 = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));        //    运算速度(因为有些循环是公共的)会有点点提高,额外的内存占用则可以忽略

    if ((B1 == NULL) || (B2 == NULL) || (B3 == NULL))
    {
        if (B1 != NULL) free(B1);
        if (B2 != NULL) free(B2);
        if (B3 != NULL) free(B3);
        return IM_STATUS_OUTOFMEMORY;
    }

    Status = IM_ExpBlur(Src, B1, Width, Height, Stride, Radius);                        
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    Status = IM_ExpBlur(Src, B2, Width, Height, Stride, Radius * 2);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    Status = IM_ExpBlur(Src, B3, Width, Height, Stride, Radius * 4);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    for (int Y = 0; Y < Height * Stride; Y++)
    {
        int DiffB1 = Src[Y] - B1[Y];
        int DiffB2 = B1[Y] - B2[Y];
        int DiffB3 = B2[Y] - B3[Y];
        Dest[Y] = IM_ClampToByte(((4 - 2 * IM_Sign(DiffB1)) * DiffB1 + 2 * DiffB2 + DiffB3) / 4 + Src[Y]);
    }

FreeMemory:
    free(B1);
    free(B2);
    free(B3);
    return Status;
}

  为避免浮点计算,我们将W1、W2、W3放大四倍,累加玩后在除以4就可以了,整个的书写过程就是按照公式(13)进行的。

其中IM_Sign函数定义如下:

inline int IM_Sign(int X)
{
    return (X >> 31) | (unsigned(-X)) >> 31;        
}

                处理前

                      处理后(半径5)

处理效果由上面两幅图的比较来说,还是相当明显的。

上面的代码中我用的ExpBlur代替了高斯模糊,关于指数模糊可以参考:SSE图像算法优化系列五:超高速指数模糊算法的实现和优化(10000*10000在100ms左右实现) 一文,他的效果和高斯模糊差不多,速度要快不少。

当指数模糊使用SSE优化后,剩下的代码用纯C实现,对于1080P的24位彩色图,测试时间为73毫秒,如果除去3次指数模糊,纯C部分代码耗时约40 ms,所以有很大的优化空间,我们就用SSE来处理。

第一,我们需要来处理IM_Sign这个函数,这个函数当参数大于0时,返回1,参数小于0时,返回-1,参数等于0时,返回0,SSE没有直接这样的函数,幸好之前收集过一个这个的自定义函数,写的也很巧妙:

inline __m128i _mm_sgn_epi16(__m128i v)
{
#ifdef __SSSE3__
    v = _mm_sign_epi16(_mm_set1_epi16(1), v); // use PSIGNW on SSSE3 and later
#else
    v = _mm_min_epi16(v, _mm_set1_epi16(1));  // use PMINSW/PMAXSW on SSE2/SSE3.
    v = _mm_max_epi16(v, _mm_set1_epi16(-1));
    //_mm_set1_epi16(1) = _mm_srli_epi16(_mm_cmpeq_epi16(v, v), 15);
    //_mm_set1_epi16(-1) = _mm_cmpeq_epi16(v, v);

#endif
    return v;
}

  如上所示,当系统只支持SSE2以及其下的版本时,直接用_mm_min_epi16和_mm_max_epi16这样的函数硬实现,这个逻辑也很好理解。

如果系统支持SSE3及其以上的版本,系统提供了_mm_sign_epi16这个函数,关于这个函数其作用解释如下:

//  extern __m128i _mm_sign_epi16 (__m128i a, __m128i b); 
//  Negate packed words in a if corresponding sign in b is less than zero. 
//  Interpreting a, b, and r as arrays of signed 16-bit integers: 
for (i = 0; i < 8; i++)
{
    if (b[i] < 0)
    {
        r[i] = -a[i];
    }
    else if (b[i] == 0)
    {
        r[i] = 0;
    }
    else
    {
        r[i] = a[i];
    }
}

  如果参数a传值为1,不就能实现IM_Sign这个效果了吗,真好玩。

解决了IM_Sign这个函数,其他的部分都很简单了,考虑到数据范围,把字节数据扩展为signed short类型处理就可以了,这样一次性可以处理8个字节的数据,修改后的代码如下所示:

int IM_MultiScaleSharpen(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                                return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                                    return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3) && (Channel != 4))                return IM_STATUS_INVALIDPARAMETER;
    int Status = IM_STATUS_OK;
    unsigned char *B1 = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));        //    如果最后的大图在这里内存有问题,一种解决办法就是下面的模糊用BoxBlur
    unsigned char *B2 = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));        //    然后不要调用现有的BoxBlur函数,而是直接在本函数实现三种不同半径的模糊
    unsigned char *B3 = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));        //    运算速度(因为有些循环是公共的)会有点点提高,额外的内存占用则可以忽略
    if ((B1 == NULL) || (B2 == NULL) || (B3 == NULL))
    {
        if (B1 != NULL) free(B1);
        if (B2 != NULL) free(B2);
        if (B3 != NULL) free(B3);
        return IM_STATUS_OUTOFMEMORY;
    }
    Status = IM_ExpBlur(Src, B1, Width, Height, Stride, Radius);                        
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    Status = IM_ExpBlur(Src, B2, Width, Height, Stride, Radius * 2);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    Status = IM_ExpBlur(Src, B3, Width, Height, Stride, Radius * 4);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    int BlockSize = 8, Block = (Height * Stride) / BlockSize;
    __m128i Zero = _mm_setzero_si128();
    __m128i Four = _mm_set1_epi16(4);
    for (int Y = 0; Y < Block * BlockSize; Y += BlockSize)
    {
        __m128i SrcV = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Src + Y)), Zero);
        __m128i SrcB1 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(B1 + Y)), Zero);
        __m128i SrcB2 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(B2 + Y)), Zero);
        __m128i SrcB3 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(B3 + Y)), Zero);
        __m128i DiffB1 = _mm_sub_epi16(SrcV, SrcB1);
        __m128i DiffB2 = _mm_sub_epi16(SrcB1, SrcB2);
        __m128i DiffB3 = _mm_sub_epi16(SrcB2, SrcB3);
        __m128i Offset = _mm_srai_epi16(_mm_add_epi16(_mm_add_epi16(_mm_mullo_epi16(_mm_sub_epi16(Four, _mm_slli_epi16(_mm_sgn_epi16(DiffB1), 1)), DiffB1), _mm_slli_epi16(DiffB2, 1)), DiffB3), 2);
        _mm_storel_epi64((__m128i *)(Dest + Y), _mm_packus_epi16(_mm_add_epi16(SrcV, Offset), Zero));
    }
    for (int Y = Block * BlockSize; Y < Height * Stride; Y++)
    {
        int DiffB1 = Src[Y] - B1[Y];
        int DiffB2 = B1[Y] - B2[Y];
        int DiffB3 = B2[Y] - B3[Y];
        Dest[Y] = IM_ClampToByte(((4 - 2 * IM_Sign(DiffB1)) * DiffB1 + 2 * DiffB2 + DiffB3) / 4 + Src[Y]);
    }
FreeMemory:
    free(B1);
    free(B2);
    free(B3);
    return Status;
}

  基本上就是按照C语言或者公式(13)所示的流程一步一步的编写,只不过把有些乘法变成了移位。

 对于1080P的彩色图像,上述改动后处理时间变为了35ms,纯C语言部分的耗时约在11ms左右,同之前的相比速度提高了4倍多,提速还是相当的明显的。

在仔细观察,觉得在IM_Sign这个的处理上还是有问题,虽然上述代码能完美解决问题,但是总觉得有改进空间,当我们把IM_Sign(DiffB1) * DiffB1放在一起观察时,就会发现这个整体不是可以直接使用_mm_sign_epi16予以实现吗,比如 _mm_sign_epi16(DiffB1, DiffB1) 难道不是吗? 这样就少了一次乘法。

最后,我们把DiffB2, DiffB3展开,合并掉一些同类项,然后把乘以相同系数的变量在合并,又可以优化掉一些计算,最终的SSE部分代码如下:

    int BlockSize = 8, Block = (Height * Stride) / BlockSize;
    __m128i Zero = _mm_setzero_si128();
    for (int Y = 0; Y < Block * BlockSize; Y += BlockSize)
    {
        __m128i SrcV = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Src + Y)), Zero);
        __m128i SrcB1 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(B1 + Y)), Zero);
        __m128i SrcB2 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(B2 + Y)), Zero);
        __m128i SrcB3 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(B3 + Y)), Zero);
        __m128i DiffB1 = _mm_sub_epi16(SrcV, SrcB1);
        __m128i Offset = _mm_add_epi16(_mm_srai_epi16(_mm_sub_epi16(_mm_slli_epi16(_mm_sub_epi16(SrcB1, _mm_sign_epi16(DiffB1, DiffB1)), 1), _mm_add_epi16(SrcB2, SrcB3)), 2), DiffB1);
        _mm_storel_epi64((__m128i *)(Dest + Y), _mm_packus_epi16(_mm_add_epi16(SrcV, Offset), Zero));
    }

  虽然测试表明,速度没有提高多少(主要是计算量太少),但这样写明显合理很多。

测试工程的地址:http://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏一心无二用,本人只专注于基础图像算法的实现与优化。

SSE图像算法优化系列七:基于SSE实现的极速的矩形核腐蚀和膨胀(最大值和最小值)算法。

  因未测试其他作者的算法时间和效率,本文不敢自称是最快的,但是速度也可以肯定说是相当快的,在一台I5机器上占用单核的资源处理 3000 * 2000的灰度...

2569
来自专栏数据小魔方

创意滑珠图!

今天要给大家分享的是一种非常有趣的滑珠图! ▽ 本文要讲解的滑珠图做法,稍微有点复杂。不过这种滑珠图在数据对比展示中,效果奇佳。小魔方参考多处教程和资料,终于还...

2544
来自专栏菩提树下的杨过

Flash/Flex学习笔记(57):实用技巧

布朗运动: varnumDots:uint=50; varfriction:Number=0.9; vardots:Array; varlife:uint=0;...

23910
来自专栏葡萄城控件技术团队

WPF/Silverlight Layout 系统概述——Measure

前言 在WPF/Silverlight当中,如果已经存在的Element无法满足你特殊的需求,你可能想自定义Element,那么就有可能会面临重写Measure...

1678
来自专栏落影的专栏

Metal入门教程(四)灰度计算

Metal入门教程(一)图片绘制 Metal入门教程(二)三维变换 Metal入门教程(三)摄像头采集渲染

1473
来自专栏郭艺帆的专栏

庖丁解牛:GIF

GIF 是一种使用 LZW 压缩,支持多张图像的容器。支持256色,透明通道为1bit。作为互联网表情包的载体,GIF 这项80年代的技术依然生生不息。

1700
来自专栏腾讯Bugly的专栏

Android自绘动画实现与优化实战——以Tencent OS录音机波形动画为实例

前言 我们所熟知的,Android 的图形绘制主要是基于 View 这个类实现。 每个 View 的绘制都需要经过 onMeasure、onLayout、onD...

3034
来自专栏申龙斌的程序人生

零基础学编程015:画些有趣的图案

从《零基础学编程014:小海龟做画》中我们学会了基本的做图命令,只需要用上循环语句,就可以画出比较复杂的图案来,比如: from turtle import *...

3049
来自专栏一心无二用,本人只专注于基础图像算法的实现与优化。

二值图像中封闭孔洞的高效填充算法(附源码)。

     鉴于心情不好,这篇文章只是简单的说说这个算法的过程。      在对图像二值化后,不管用的是什么二值算法,总会存在一些瑕疵,这个时候我们就需要进行一些...

3096
来自专栏生信宝典

R语言学习 - 热图美化

热图美化 上一期的绘图命令中,最后一行的操作抹去了之前设定的横轴标记的旋转,最后出来的图比较难看。 上次我们是这么写的 p <- p + xlab("samp...

2998

扫码关注云+社区