2015年龚博士的曲率滤波算法刚出来的时候，在图像处理界也曾引起不小的轰动，特别是其所说的算法的简洁性，以及算法的效果、执行效率等方面较其他算法均有一定的优势，我在该算法刚出来时也曾经有关注，不过那个时候看到是迭代的算法，而且迭代的次数还蛮多了，就觉得算法应该不会太快，所以就放弃了对其进一步优化。最近，又偶尔一次碰触到该文章和代码，感觉还是有蛮大的优化空间的，所以抽空简单的实现他的算法。 　 该算法作者已经完全开源，项目地址见：https://github.com/YuanhaoGong/CurvatureFilter ， 里面有matlab\c++\python等语言的代码，其中matlab的代码比较简洁，C++的是基于opencv的，而且里面包含了很多其他方面的代码，整体看起来我感觉有点乱。我只是稍微浏览了下。 　　通读matlab的代码，其和论文里提供的伪代码好像除了TV算法外，其他的都基本对应不上，我不知道是怎么回事，因此，本文仅以TV算法的优化作为例子，而且TV算法在GC\MC\TV算法中属于实现最为复杂和耗时的一个，也是最能反映优化极限的例子， 因此处理好了这个算法，另外两个的优化就是水道渠成的事情了。 　　算法的原理我不介绍，我也不在行，论文中给出的伪代码如下所示： 　　 　　简单的说，对于每个像素，我们以其为中心，领域3*3的点，做算法15的操作，新的像素值就是原始像素值加上dm。但是我们做的时候是把整副图像分成四块，如上图右侧图所示，分为白色圆、黑色圆、白色三角形、黑色三角形，这四块独立更新，更新完后的心像素值可以为尚未更新的像素所使用，一副图全部更新后，再次迭代该过程，直到达到指定迭代次数为止。 　 在国内网站上搜索，我发现网友CeleryChen在他的博文中曲率滤波的简单实现只有源码提到了该算法的实现代码，我看了下，和我的书写习惯还是很类似的（比原作者的代码看起来舒服多了），不过他写的不是TV算法， 我修改为TV算法后，核心代码如下所示： 复制代码 void TV_Curvature_Filter(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int StartRow, int StartCol) { for (int Y = StartRow; Y < Height - 1; Y += 2) { unsigned char *LinePC = Src + Y * Stride; unsigned char *LinePL = Src + (Y - 1) * Stride; unsigned char *LinePN = Src + (Y + 1) * Stride; unsigned char *LinePD = Dest + Y * Stride; for (int X = StartCol; X < Width - 1; X += 2) { short Dist[8] = { 0 }; int C5 = 5 * LinePC[X]; Dist[0] = (LinePL[X - 1] + LinePL[X] + LinePC[X - 1] + LinePN[X - 1] + LinePN[X]) - C5; Dist[1] = (LinePL[X] + LinePL[X + 1] + LinePC[X + 1] + LinePN[X] + LinePN[X + 1]) - C5; Dist[2] = (LinePL[X - 1] + LinePL[X] + LinePL[X + 1] + LinePC[X - 1] + LinePC[X + 1]) - C5; Dist[3] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePC[X - 1] + LinePC[X + 1]) - C5; Dist[4] = (LinePL[X - 1] + LinePL[X] + LinePL[X + 1] + LinePC[X - 1] + LinePN[X - 1]) - C5; Dist[5] = (LinePL[X - 1] + LinePL[X + 1] + LinePL[X + 1] + LinePC[X + 1] + LinePN[X + 1]) - C5; Dist[6] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePL[X - 1] + LinePC[X - 1]) - C5; Dist[7] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePL[X + 1] + LinePC[X + 1]) - C5; __m128i Dist128 = _mm_loadu_si128((const __m128i*)Dist); __m128i AbsDist128 = _mm_abs_epi16(Dist128); __m128i MinPos128 = _mm_minpos_epu16(AbsDist128); LinePD[X] = IM_ClampToByte(LinePC[X] + Dist[_mm_extract_epi16(MinPos128, 1)] / 5) ; } } } void IM_CurvatureFilter(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Iteration) { unsigned char* Temp = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char)); memcpy(Temp, Src, Height * Stride * sizeof(unsigned char)); memcpy(Dest, Src, Height * Stride * sizeof(unsigned char)); for (int Y = 0; Y < Iteration; Y++) { TV_Curvature_Filter(Temp, Dest, Width, Height, Stride, 1, 1); memcpy(Temp, Dest, Height * Stride * sizeof(unsigned char)); TV_Curvature_Filter(Temp, Dest, Width, Height, Stride, 2, 2); memcpy(Temp, Dest, Height * Stride * sizeof(unsigned char)); TV_Curvature_Filter(Temp, Dest, Width, Height, Stride, 1, 2); memcpy(Temp, Dest, Height * Stride * sizeof(unsigned char)); TV_Curvature_Filter(Temp, Dest, Width, Height, Stride, 2, 1); memcpy(Temp, Dest, Height * Stride * sizeof(unsigned char)); } free(Temp); } 复制代码 　　代码非常简单和简短，确实不超过100行，我们将迭代次数设置为50次，对比对应的matlab效果如下： 　　　　　　　　　　　　　　　　　　原图　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 　 　　C++代码的效果（迭代50次）　　　　　　　　　　　　　　　　　　　　 　　　　matlab代码的效果 　　仔细观察，发现C++代码的效果在头发，帽子部位明显和matlab的不同，所以难怪CeleryChen在代码注释里有这样一段话：打算自己实现一个更快更好的曲率滤波的算法将其开源, 但发现按照作者的博士论文的算法出来的效果,达不到作者文献中的效果。 　　其实不是CeleryChen的代码逻辑问题，而是他可能没有注意到数据范围，在上述代码中，他使用unsigned char作为中间数据，而在更新的累加过程中使用了类似这样的语句： 　　　　　　LinePD[X] = IM_ClampToByte(LinePC[X] + Dist[_mm_extract_epi16(MinPos128, 1)]); 　　Clamp函数把数据再次抑制到0和255之间，这个在单次处理时问题不大，但是在迭代过程中就会出现较大的问题，总体来说同样的迭代次数效果要明显弱于matlab的代码，所以只要这个数据类型进行扩展，就可以，比如更改为浮点类型，代码如下所示： 复制代码 void TV_Curvature_Filter(float *Src, float *Dest, int Width, int Height, int Stride, int StartRow, int StartCol) { for (int Y = StartRow; Y < Height - 1; Y += 2) { float *LinePC = Src + Y * Stride; float *LinePL = Src + (Y - 1) * Stride; float *LinePN = Src + (Y + 1) * Stride; float *LinePD = Dest + Y * Stride; for (int X = StartCol; X < Width - 1; X += 2) { float Dist[8] = { 0 }; float C5 = 5 * LinePC[X]; Dist[0] = (LinePL[X - 1] + LinePL[X] + LinePC[X - 1] + LinePN[X - 1] + LinePN[X]) - C5; Dist[1] = (LinePL[X] + LinePL[X + 1] + LinePC[X + 1] + LinePN[X] + LinePN[X + 1]) - C5; Dist[2] = (LinePL[X - 1] + LinePL[X] + LinePL[X + 1] + LinePC[X - 1] + LinePC[X + 1]) - C5; Dist[3] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePC[X - 1] + LinePC[X + 1]) - C5; Dist[4] = (LinePL[X - 1] + LinePL[X] + LinePL[X + 1] + LinePC[X - 1] + LinePN[X - 1]) - C5; Dist[5] = (LinePL[X - 1] + LinePL[X + 1] + LinePL[X + 1] + LinePC[X + 1] + LinePN[X + 1]) - C5; Dist[6] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePL[X - 1] + LinePC[X - 1]) - C5; Dist[7] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePL[X + 1] + LinePC[X + 1]) - C5; float Min = abs(Dist[0]); int Index = 0; for (int Z = 1; Z < 8; Z++) { if (abs(Dist[Z]) < Min) { Min = abs(Dist[Z]); Index = Z; } } LinePD[X] = LinePC[X] + Dist[Index] * 0.2f; } } } void IM_CurvatureFilter(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Iteration) { float *Temp1 = (float *)malloc(Height * Stride * sizeof(float)); float *Temp2 = (float *)malloc(Height * Stride * sizeof(float)); for (int Y = 0; Y < Height * Stride; Y++) Temp1[Y] = Src[Y]; memcpy(Temp2, Temp1, Height * Stride * sizeof(float)); for (int Y = 0; Y < Iteration; Y++) { TV_Curvature_Filter(Temp1, Temp2, Width, Height, Stride, 1, 1); memcpy(Temp1, Temp2, Height * Stride * sizeof(float)); TV_Curvature_Filter(Temp1, Temp2, Width, Height, Stride, 2, 2); memcpy(Temp1, Temp2, Height * Stride * sizeof(float)); TV_Curvature_Filter(Temp1, Temp2, Width, Height, Stride, 1, 2); memcpy(Temp1, Temp2, Height * Stride * sizeof(float)); TV_Curvature_Filter(Temp1, Temp2, Width, Height, Stride, 2, 1); memcpy(Temp1, Temp2, Height * Stride * sizeof(float)); } for (int Y = 0; Y < Height * Stride; Y++) Dest[Y] = IM_ClampToByte((int)(Temp2[Y] + 0.4999999f)); free(Temp1); free(Temp2); } 复制代码 　　结果如下： 　　　　　　　　　　 　　C++代码效果　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 matlab代码效果　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 　　 差异图 　　差异图不是全黑，这个和m代码里局部细节和C++有点不同有关。 　　到此，我们彻底消除CeleryChen所说的的算法效果问题。 　　论文中说到了讲图像分成四个块，如上述的代码所示，分四个部分更新，那其实这里就有个问题，四个块的更新顺序是什么，或者说如果你指定一个顺序，那你的理论依据是啥，因为一个块更新后，后面一个块的更新会依赖于前面已经更新的块的数据，所以这个顺序对结果可定是有影响的，我在论文里没有找到相关说法，那么好，我就实践，我把不同可能的更新顺序的结果都做出来， 结果发现，大家确实有区别，但是区别都是在一两个像素之间，因此，这个顺序就显得不是很重要了。 　　但是，我们不分块，又会存在什么问题呢，理论上有没有问题，我不知道，但是我们就实测下，不分块，代码如下所示： 复制代码 void TV_Curvature_Filter(float *Src, float *Dest, int Width, int Height, int Stride) { for (int Y = 1; Y < Height - 1; Y++) { float *LinePC = Src + Y * Stride; float *LinePL = Src + (Y - 1) * Stride; float *LinePN = Src + (Y + 1) * Stride; float *LinePD = Dest + Y * Stride; for (int X = 1; X < Width - 1; X++) { float Dist[8] = { 0 }; float C5 = 5 * LinePC[X]; Dist[0] = (LinePL[X - 1] + LinePL[X] + LinePC[X - 1] + LinePN[X - 1] + LinePN[X]) - C5; Dist[1] = (LinePL[X] + LinePL[X + 1] + LinePC[X + 1] + LinePN[X] + LinePN[X + 1]) - C5; Dist[2] = (LinePL[X - 1] + LinePL[X] + LinePL[X + 1] + LinePC[X - 1] + LinePC[X + 1]) - C5; Dist[3] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePC[X - 1] + LinePC[X + 1]) - C5; Dist[4] = (LinePL[X - 1] + LinePL[X] + LinePL[X + 1] + LinePC[X - 1] + LinePN[X - 1]) - C5; Dist[5] = (LinePL[X - 1] + LinePL[X + 1] + LinePL[X + 1] + LinePC[X + 1] + LinePN[X + 1]) - C5; Dist[6] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePL[X - 1] + LinePC[X - 1]) - C5; Dist[7] = (LinePN[X - 1] + LinePN[X] + LinePN[X + 1] + LinePL[X + 1] + LinePC[X + 1]) - C5; float Min = abs(Dist[0]); int Index = 0; for (int Z = 1; Z < 8; Z++) { if (abs(Dist[Z]) < Min) { Min = abs(Dist[Z]); Index = Z; } } LinePD[X] = LinePC[X] + Dist[Index] * 0.2f; } } } void IM_CurvatureFilter_01(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Iteration) { float *Temp1 = (float *)malloc(Height * Stride * sizeof(float)); float *Temp2 = (float *)malloc(Height * Stride * sizeof(float)); for (int Y = 0; Y < Height * Stride; Y++) Temp1[Y] = Src[Y]; memcpy(Temp2, Temp1, Height * Stride * sizeof(float)); for (int Y = 0; Y < Iteration; Y++) { TV_Curvature_Filter(Temp1, Temp2, Width, Height, Stride); memcpy(Temp1, Temp2, Height * Stride * sizeof(float)); } for (int Y = 0; Y < Height * Stride; Y++) Dest[Y] = IM_ClampToByte((int)(Temp2[Y] + 0.4999999f)); free(Temp1); free(Temp2); } 复制代码 　　注意到在TV_Curvature_Filter函数中的两重for循环里每次循环变量自增值已经修改了，如果分块则一次增加2个像素，同时看到在前面未分块的代码汇中循环结束的标记是宽度和高度都减一，这个其实只是为了编码方便，放置访问越位内存，在matlab代码中，作者是做了扩边处理的，这些都不重要。 　　我们来看看处理结果的比较： 　 　　　　　　　　　 　不分块结果图　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 　分块结果图　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　差异图 　　确实有些差异，并且差异的比各分块之间的差异要大， 但是我觉得没有到那种已经有质的区别的档次。因此后续的优化我是基于不分块的来描述了，但是这种优化也是可以稍作修改就应用于分块的优化（代码会显得稍微繁琐一点，速度也会稍有下降）。 　　我们首先记录下目前的速度，不分块50次迭代512*512的标准lena灰度图耗时约 390ms（VS2013编译器，默认启用了SSE优化）。 　　第一个小优化：每次迭代里有个memcpy函数，目的是把中间的临时目的数据拷贝到临时源数据中，这里其实可以不用这个拷贝，直接使用个指针交换就可以了，如下所示： 复制代码 void IM_CurvatureFilter_01(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Iteration) { float *Temp1 = (float *)malloc(Height * Stride * sizeof(float)); float *Temp2 = (float *)malloc(Height * Stride * sizeof(float)); float *Temp = NULL; for (int Y = 0; Y < Height * Stride; Y++) Temp1[Y] = Src[Y]; memcpy(Temp2, Temp1, Height * Stride * sizeof(float)); for (int Y = 0; Y < Iteration; Y++) { TV_Curvature_Filter(Temp1, Temp2, Width, Height, Stride); Temp = Temp1, Temp1 = Temp2, Temp2 = Temp; //memcpy(Temp1, Temp2, Height * Stride * sizeof(float)); } for (int Y = 0; Y < Height * Stride; Y++) Dest[Y] = IM_ClampToByte((int)(Temp2[Y] + 0.4999999f)); free(Temp1); free(Temp2); } 复制代码 　　执行速度变为约385ms，这只是个小菜，没啥意思。 　　第二步尝试，改变下算法的逻辑，我们注意到在Dist的8个元素计算中，其实是有很多重复计算的，我们如果不管Dist计算的顺序，其实是可以总结出其运算规律： 　　如果我们计算出第一个Dist中五个元素的和，后续就只要加一个进入的元素并且减去一个离开的元素就可以了，而且每次的减C5也可以只在第一个里做，这样TV_Curvature_Filter变为如下形式： 复制代码 void TV_Curvature_Filter(float *Src, float *Dest, int Width, int Height, int Stride) { for (int Y = 1; Y < Height - 1; Y++) { float *LinePC = Src + Y * Stride; float *LinePL = Src + (Y - 1) * Stride; float *LinePN = Src + (Y + 1) * Stride; float *LinePD = Dest + Y * Stride; for (int X = 1; X < Width - 1; X++) { float Dist[8] = { 0 }; Dist[0] = (LinePL[X - 1] + LinePL[X] + LinePL[X + 1] + LinePC[X + 1] + LinePN[X + 1]) - 5 * LinePC[X]; Dist[1] = Dist[0] - LinePL[X - 1] + LinePN[X]; Dist[2] = Dist[1] - LinePL[X] + LinePN[X - 1]; Dist[3] = Dist[2] - LinePL[X + 1] + LinePC[X - 1]; Dist[4] = Dist[3] - LinePC[X + 1] + LinePL[X - 1]; Dist[5] = Dist[4] - LinePN[X + 1] + LinePL[X]; Dist[6] = Dist[5] - LinePN[X] + LinePL[X + 1]; Dist[7] = Dist[6] - LinePN[X - 1] + LinePC[X + 1]; float Min = abs(Dist[0]); int Index = 0; for (int Z = 1; Z < 8; Z++) { if (abs(Dist[Z]) < Min) { Min = abs(Dist[Z]); Index = Z; } } LinePD[X] = LinePC[X] + Dist[Index] * 0.2f; } } } 复制代码 　　核心的计算由原先的32次加法及8次减法变化为11次加法及8次减法，计算大为减少，这应该能加速不少吧。 　　按下F5，我靠运行时间变为了420ms，真是大失所望，为什么计算量少了但运行时间却多了呢，反汇编看编译器生成的代码，确实是修改的指令多很多，修改后少了不少。我自己想了想，觉得主要还是因为修改后的算法的计算是前后依赖的，在没算法前面的指令后续的指令是无法执行的，而编译器的在指令级别也是可以实现指令级并行，修改后的方案就无法重复利用这个特性了。 　　因此，这条路实际上是个失败的方案，龚博士在matlab代码里也有类似的优化方式，我把他修改为C代码后，得到的结论也是一样的。但是他 并不是一无是处，后面还是说到。 　　还有一些常规的优化，比如内部的8次循环展开，不用数组，直接用变量代替中间的距离变量等等，实践都表面没有啥效果。因此，我们不浪费时间，直接进行SSE优化了。 　　针对上述代码进行SSE自行优化应该说不是很困难的事情，可以看到，各个点的计算是独立的，和前后像素之间的计算是毫无干扰的，这种情况非常适合并行化（不管是SSE还是GPU都是一样的道理），基本上只要把对应的普通C运算符转换对一个的SIMD指令就OK了，那我们稍微关注下几个SIMD没有的运算符。 　　第一个是abs， SSE提供了大量的求绝对值的指令，但唯独对浮点数没有，看不懂这个世界，不过也没关系，浮点数普通C语言我们有个常用的方式实现绝对值计算如下所示： 复制代码 inline float IM_AbsF(float V) { int I = ((*(int *)&V) & 0x7fffffff); return (*(float *)&I); } 复制代码 　　有这个，对应的__m128数据类型计算也很简单： 复制代码 // 浮点数据的绝对值计算 inline __m128 _mm_abs_ps(__m128 v) { static const int i = 0x7fffffff; float mask = *(float*)&i; return