超越halcon速度的二值图像的腐蚀和膨胀，实现目前最快的半径相关类算法（附核心源码）。
我在两年前的博客里曾经写过 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菜单。
