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

时间:2020-12-09 14:13:19

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

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

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

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

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

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

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

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

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

 int BlockSize = , Block = (Width - ) / BlockSize;
__m128i OldSum = _mm_set1_epi32(LastSum);
__m128 Inv128 = _mm_set1_ps(Inv);
for (int X = ; X < Block * BlockSize + ; X += BlockSize)
{
__m128i ColValueOut = _mm_loadu_si128((__m128i *)(ColValue + X - ));
__m128i ColValueIn = _mm_loadu_si128((__m128i *)(ColValue + X + Radius + Radius));
__m128i ColValueDiff = _mm_sub_epi32(ColValueIn, ColValueOut); // P3 P2 P1 P0
__m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, )); // P3+P2 P2+P1 P1+P0 P0
__m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, )); // P3+P2+P1+P0 P2+P1+P0 P1+P0 P0
__m128i NewSum = _mm_add_epi32(OldSum, Value);
OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(, , , )); // 重新赋值为最新值
__m128 Mean = _mm_mul_ps(_mm_cvtepi32_ps(NewSum), Inv128);
_mm_storesi128_4char(_mm_cvtps_epi32(Mean), LinePD + X);
}

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

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

LinePD[X + ] = NewSum.m128i_i32[] >  ?  : ;
LinePD[X + ] = NewSum.m128i_i32[] > ? : ;
LinePD[X + ] = NewSum.m128i_i32[] > ? : ;
LinePD[X + ] = NewSum.m128i_i32[] > ? : ;

  以上代码只是示意,如果真的这样写,会破坏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, , Width * sizeof(unsigned char));
for (int Z = -Radius; Z <= Radius; Z++)
{
unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride;
int BlockSize = , Block = Width / BlockSize;
for (int X = ; 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] == ? : ); // 更新列数据
}
}

  普通的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菜单。

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