超越halcon速度的二值图像的腐蚀和膨胀，实现目前最快的半径相关类算法（附核心源码）

2021-07-11 19:59:18 阅读：229 来源： 互联网

标签：__ mm halcon 算法源码 m128i SSE 二值 255

本篇博文来自博主Imageshop，打赏或想要查阅更多内容可以移步至Imageshop。

转载自：https://www.cnblogs.com/Imageshop/p/10563354.html 侵删

　我在两年前的博客里曾经写过 SSE图像算法优化系列七：基于SSE实现的极速的矩形核腐蚀和膨胀（最大值和最小值）算法一文，通过SSE的优化把矩形核心的腐蚀和膨胀做到了不仅和半径无关，而且速度也相当的快，当时在被博文的评论里有博友提出了如下的问题：

#1楼
2018-02-21 20:26 | 胡一谭  
博主的思路很巧妙，只是这个算法本身还是不够快，优化效果与商业软件还是有比较大差距，4096X8192大小的的灰度图商业软件（halcon）只需要33ms， 本文需要250ms，考虑到商业软件采用多核优化，我测试机器是4核，
通常优化加速比在3倍左右，因此，本文并行化后的理论耗时为250/3=83.33ms。但我采用OpenMP对本文算法进行优化后达不到3倍的加速比。还是需要寻找更好的思路。

　　当时看到这个评论后，真的觉得这博友是不是搞错了，这么大的图像，怎么可能只要33ms就处理完了呢，就是最简单的一个图像处理算法，反色（Invert）经过极度优化后也需要大概7/8毫秒的，所以我当时内心是不认可这个速度的。

　　后续我也在考虑二值图像的这个特殊性，曾经有考虑过比如膨胀时，遇到有个是白色的像素则停止循环，也考虑过使用直方图的方式进行优化，毕竟直方图也只有两个像素了，但是也还是达不到上述速度，有些甚至还更慢。所以后续一直也没有什么进步。

　　前几日，网友LQC-Jack突然又再次提到了这个问题，他认为针对这个问题确实有更快的方法，毕竟二值得特殊性摆在那里：

　　其中的“你box滤波的，sum>0当前点就是255” 这个是关键，是啊，针对二值图求局部矩形内的最大值，和求二值图像的局部均值如果我们能够建立起联系，那么就可以借助于快速的局部均值算法间接的实现腐蚀或膨胀，我在博客里有多篇文章提到了局部均值的终极优化，特别是SSE图像算法优化系列十三：超高速BoxBlur算法的实现和优化（Opencv的速度的五倍）一文中提到的方式，效率及其高，针对4096X8192的二值图也就是30ms左右能搞定。希望燃起。

　　那如何将两者搭桥呢，仔细想想确实很简单，如果是求最大值（膨胀），那么只要局部有一个像素为255，结果就为255，此时的局部均值必然大于0 （考虑实际因素，应该是局部累加值，因为考虑最后的整除，不排除某个局部区域，只有一个白点，当局部过大时，整除后的结果可能也为0），而只有所有局部内的像素都为0是，最大值才为0，这个时候的局部累加值也必然为0。如果是求最,小值（腐蚀），只要局部有一个像素为0值，结果就为0，只有局部所有像素都为255，结果才为255，那么这里的信息反馈到局部均值就等同于说平均值为255，则结果为255，否则结果就为0（同样的道理，这里实际编程时要用局部累加值，而不是平均值）。

　　如此一来，我们会发现，这种实现过程相比标准的方框模糊来说还少了一些步骤，我们先贴下我SSE优化方框模糊的核心部分：

复制代码

 1 int BlockSize = 4, Block = (Width - 1) / BlockSize;
 2 __m128i OldSum = _mm_set1_epi32(LastSum);
 3 __m128 Inv128 = _mm_set1_ps(Inv);
 4 for (int X = 1; X < Block * BlockSize + 1; X += BlockSize)
 5 {
 6     __m128i ColValueOut = _mm_loadu_si128((__m128i *)(ColValue + X - 1));
 7     __m128i ColValueIn = _mm_loadu_si128((__m128i *)(ColValue + X + Radius + Radius));
 8     __m128i ColValueDiff = _mm_sub_epi32(ColValueIn, ColValueOut);                            //    P3 P2 P1 P0                                                
 9     __m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4));        //    P3+P2 P2+P1 P1+P0 P0
10     __m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8));                 //    P3+P2+P1+P0 P2+P1+P0 P1+P0 P0
11     __m128i NewSum = _mm_add_epi32(OldSum, Value);
12     OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3));                              //    重新赋值为最新值
13     __m128 Mean = _mm_mul_ps(_mm_cvtepi32_ps(NewSum), Inv128);
14     _mm_storesi128_4char(_mm_cvtps_epi32(Mean), LinePD + X);
15 }

复制代码

　　注意第14及第15行为求均值并最终保存数据到内存中的过程，其中的NewSum中保存即为累加的值，注意这里因为有除法，所以借用了浮点版本的相关指令，同时增加了相关的类型转换过程。

　　这里，为了适应我们的腐蚀和膨胀的需求，这两句是不需要的，按照上述分析，比如膨胀效果，只需作如下改动：

LinePD[X + 0] = NewSum.m128i_i32[0] > 0 ? 255 : 0;
LinePD[X + 1] = NewSum.m128i_i32[1] > 0 ? 255 : 0;
LinePD[X + 2] = NewSum.m128i_i32[2] > 0 ? 255 : 0;
LinePD[X + 3] = NewSum.m128i_i32[3] > 0 ? 255 : 0;

　　以上代码只是示意，如果真的这样写，会破坏SSE算法的整体的和谐性，而且这种SSE中穿插普通C代码会带来性能上的极大损失，一种处理方法如下所示：

__m128i Flag = _mm_cmpgt_epi32(NewSum, _mm_setzero_si128());
Flag = _mm_packs_epi32(Flag, Flag);
*((int *)(LinePD + X)) = _mm_cvtsi128_si32(_mm_packs_epi16(Flag, Flag));

　　我们利用SSE中比较运算符的特殊性，产生诸如0XFFFFFFFF这样的结果，然后在通过有关饱和运算将他们减低到8位，注意上面使用的都是有符号的饱和计算。

　　对于腐蚀的过程，你知道怎么写吗？

　　我们经过简单测试，处理一副4096X8192大小的二值图，任意的半径大小，耗时基本稳定在24ms左右，比boxblur也快了很多。

　　我也构思过不实用累加和的方式判断，比如使用或运算或者与运算，但是都是解决不了进出像素的处理问题，因此，整体看来是还是用累加最为科学。

　　其实对于半径比较小时，还是有更为快速的方法的，这里稍微简单描述下，但是可能很多人看不懂。

　　在我们上述的实现中，我们用的是int类型的数据来保存累加值，这是因为半径稍微大一点累加值就可能超过short类型所能表达的范围，但是int类型SSE一次只能处理4个，而short类型数据SSE一次能处理8个，因此，如果做适当的变动是否有可能使用short类型呢，是用可能的。

　　因为是二值图，所以就只有0和255两个值，0值无所谓，那如果我们把255这个值修改成1，那么在半径不大于某个数值（64还是其他数，可以自己画一画）时，累加值将可控在short类型所能表达的范围。

　　这是还有个问题就是，255这个值如何变为1，如果使用_mm_blendv_epi8集合有关判断语句是可以实现的，但是这个Blend是比较耗时的，反而得不偿失。一个最好的办法就是充分利用无符号和有符号数之间的特点，当我们把一个等于255的unsigned char数据类型强制转换为signed char时，他的值就等于-1，和我们要的值1相反，这个时候我们原本代码里的_mm_add_epi8接可以使用_mm_sub_epi8代替，反之亦然。而在SSE里，这种类型转换还不需要强制进行，因为他直接操作内存。

　　我们贴下下面的代码可能有人就能明白是什么意思了。

复制代码

memset(ColValue + Radius, 0, Width * sizeof(unsigned char));
for (int Z = -Radius; Z <= Radius; Z++)
{
    unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride;
    int BlockSize = 16, Block = Width / BlockSize;
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        unsigned char *DestP = ColValue + X + Radius;
        __m128i Sample = _mm_loadu_si128((__m128i *)(LinePS + X));                //    255成为-1
        _mm_storeu_si128((__m128i *)DestP, _mm_sub_epi8(_mm_loadu_si128((__m128i *)DestP), Sample));
    }
    for (int X = Block * BlockSize; X < Width; X++)
    {
        ColValue[X + Radius] += (LinePS[X] == 255 ? 1 : 0);                                            //    更新列数据
    }
}

复制代码

　　普通的C代码部分及时直接实现，而SSE部分，并没有看到明显的255到1之间的转换，一起都在那几句简简单单的代码中。

　　通过这种相关的优化，大概4096X8192的图能做到12到13毫秒之间，已经完全超过了Halcon的速度。

　　halcon中的腐蚀和膨胀也有圆形半径的，同样的半径下圆形半径在halcon中的耗时大概是矩形半径的8倍左右，我相信halcon的圆形半径的算法也是通过EDM算法来实现的，详见SSE图像算法优化系列二十五:二值图像的Euclidean distance map（EDM)特征图计算及其优化一文，而我这里也差不都是这样的时间比例。

　　极度优化版本工程：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar，见Binary->Processing->Erode/Dilate菜单。

标签：__,mm,halcon,算法,源码,m128i,SSE,二值,255
来源： https://blog.csdn.net/weixin_39450742/article/details/118612063

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

超越halcon速度的二值图像的腐蚀和膨胀，实现目前最快的半径相关类算法（附核心源码）