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

来源：互联网发布：九次方大数据公司老板编辑：程序博客网时间：2024/05/21 07:51

无意中浏览一篇文章，中间提到了基于多尺度的图像的细节提升算法，尝试了一下，还是有一定的效果的，结合最近一直研究的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));    }

阅读全文

0 0