【算法随记三】小半径中值模糊的急速实现（16MB图7.5ms实现） + Photoshop中蒙尘和划痕算法解读。
在本人的博客里，分享了有关中值模糊的O(1)算法，详见：任意半径中值滤波（扩展至百分比滤波器）O(1)时间复杂度算法的原理、实现及效果 ，这里的算法的执行时间和参数是无关的。整体来说，虽然速度也很快，但是在某些特殊情况下我们还是需要更快的速度。特别是对于小半径的中值，我们有理由去对其进一步的优化的。本文我们进一步探讨这个问题。
一、3*3中值模糊
首先我们来看看半径为1的中值，此时涉及到的领域为3*3，共9个像素，那么最传统的实现方式就是对9个像素直接进行排序，这里我们直接使用系统的排序函数qsort，一种简单的代码如下所示：
int __cdecl ComparisonFunction(const void *X, const void *Y) // 一定要用__cdecl这个标识符
{
unsigned char Dx = *(unsigned char *)X;
unsigned char Dy = *(unsigned char *)Y;
if (Dx < Dy)
return -1;
else if (Dx > Dy)
return +1;
else
return 0;
}
void MedianBlur3X3_Ori(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
int Channel = Stride / Width;
if (Channel == 1)
{
unsigned char Array[9];
for (int Y = 1; Y < Height - 1; Y++)
{
unsigned char *LineP0 = Src + (Y - 1) * Stride + 1;
unsigned char *LineP1 = LineP0 + Stride;
unsigned char *LineP2 = LineP1 + Stride;
unsigned char *LinePD = Dest + Y * Stride + 1;
for (int X = 1; X < Width - 1; X++)
{
Array[0] = LineP0[X - 1]; Array[1] = LineP0[X]; Array[2] = LineP0[X + 1];
Array[3] = LineP1[X - 1]; Array[4] = LineP1[X]; Array[5] = LineP2[X + 1];
Array[6] = LineP2[X - 1]; Array[7] = LineP2[X]; Array[8] = LineP2[X + 1];
qsort(Array, 9, sizeof(unsigned char), &ComparisonFunction);
LinePD[X] = Array[4];
}
}
}
else
{
unsigned char ArrayB[9], ArrayG[9], ArrayR[9];
for (int Y = 1; Y < Height - 1; Y++)
{
unsigned char *LineP0 = Src + (Y - 1) * Stride + 3;
unsigned char *LineP1 = LineP0 + Stride;
unsigned char *LineP2 = LineP1 + Stride;
unsigned char *LinePD = Dest + Y * Stride + 3;
for (int X = 1; X < Width - 1; X++)
{
ArrayB[0] = LineP0[-3]; ArrayG[0] = LineP0[-2]; ArrayR[0] = LineP0[-1];
ArrayB[1] = LineP0[0]; ArrayG[1] = LineP0[1]; ArrayR[1] = LineP0[2];
ArrayB[2] = LineP0[3]; ArrayG[2] = LineP0[4]; ArrayR[2] = LineP0[5];
ArrayB[3] = LineP1[-3]; ArrayG[3] = LineP1[-2]; ArrayR[3] = LineP1[-1];
ArrayB[4] = LineP1[0]; ArrayG[4] = LineP1[1]; ArrayR[4] = LineP1[2];
ArrayB[5] = LineP1[3]; ArrayG[5] = LineP1[4]; ArrayR[5] = LineP1[5];
ArrayB[6] = LineP2[-3]; ArrayG[6] = LineP2[-2]; ArrayR[6] = LineP2[-1];
ArrayB[7] = LineP2[0]; ArrayG[7] = LineP2[1]; ArrayR[7] = LineP2[2];
ArrayB[8] = LineP2[3]; ArrayG[8] = LineP2[4]; ArrayR[8] = LineP2[5];
qsort(ArrayB, 9, sizeof(unsigned char), &ComparisonFunction);
qsort(ArrayG, 9, sizeof(unsigned char), &ComparisonFunction);
qsort(ArrayR, 9, sizeof(unsigned char), &ComparisonFunction);
LinePD[0] = ArrayB[4];
LinePD[1] = ArrayG[4];
LinePD[2] = ArrayR[4];
LineP0 += 3;
LineP1 += 3;
LineP2 += 3;
LinePD += 3;
}
}
}
}
代码很简洁和清晰，我们没有处理边缘的那一圈像素，这无关精要，我们编译后测试，结果如下所示：
1920*1080大小的24位图像，平均用时1280ms，灰度图像平均用时460ms，这相当的慢，无法接受。
下面我们稍微对其进行下提速。
对于9个数据的排序，我们可以对其进行特殊的处理，因为数据的个数是确定的，按照理论分析，没有必要进行大规模的比较，实际只需要进行19次比较就可以了。修改后算法如下所示：
inline void Swap(int &X, int &Y)
{
X ^= Y;
Y ^= X;
X ^= Y;
}
void MedianBlur3X3_Faster(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
int Channel = Stride / Width;
if (Channel == 1)
{
for (int Y = 1; Y < Height - 1; Y++)
{
unsigned char *LineP0 = Src + (Y - 1) * Stride + 1;
unsigned char *LineP1 = LineP0 + Stride;
unsigned char *LineP2 = LineP1 + Stride;
unsigned char *LinePD = Dest + Y * Stride + 1;
for (int X = 1; X < Width - 1; X++)
{
int Gray0, Gray1, Gray2, Gray3, Gray4, Gray5, Gray6, Gray7, Gray8;
Gray0 = LineP0[X - 1]; Gray1 = LineP0[X]; Gray2 = LineP0[X + 1];
Gray3 = LineP1[X - 1]; Gray4 = LineP1[X]; Gray5 = LineP2[X + 1];
Gray6 = LineP2[X - 1]; Gray7 = LineP2[X]; Gray8 = LineP2[X + 1];
if (Gray1 > Gray2) Swap(Gray1, Gray2);
if (Gray4 > Gray5) Swap(Gray4, Gray5);
if (Gray7 > Gray8) Swap(Gray7, Gray8);
if (Gray0 > Gray1) Swap(Gray0, Gray1);
if (Gray3 > Gray4) Swap(Gray3, Gray4);
if (Gray6 > Gray7) Swap(Gray6, Gray7);
if (Gray1 > Gray2) Swap(Gray1, Gray2);
if (Gray4 > Gray5) Swap(Gray4, Gray5);
if (Gray7 > Gray8) Swap(Gray7, Gray8);
if (Gray0 > Gray3) Swap(Gray0, Gray3);
if (Gray5 > Gray8) Swap(Gray5, Gray8);
if (Gray4 > Gray7) Swap(Gray4, Gray7);
if (Gray3 > Gray6) Swap(Gray3, Gray6);
if (Gray1 > Gray4) Swap(Gray1, Gray4);
if (Gray2 > Gray5) Swap(Gray2, Gray5);
if (Gray4 > Gray7) Swap(Gray4, Gray7);
if (Gray4 > Gray2) Swap(Gray4, Gray2);
if (Gray6 > Gray4) Swap(Gray6, Gray4);
if (Gray4 > Gray2) Swap(Gray4, Gray2);
LinePD[X] = Gray4;
}
}
}
else
{
for (int Y = 1; Y < Height - 1; Y++)
{
unsigned char *LineP0 = Src + (Y - 1) * Stride + 3;
unsigned char *LineP1 = LineP0 + Stride;
unsigned char *LineP2 = LineP1 + Stride;
unsigned char *LinePD = Dest + Y * Stride + 3;
for (int X = 1; X < Width - 1; X++)
{
int Blue0, Blue1, Blue2, Blue3, Blue4, Blue5, Blue6, Blue7, Blue8;
int Green0, Green1, Green2, Green3, Green4, Green5, Green6, Green7, Green8;
int Red0, Red1, Red2, Red3, Red4, Red5, Red6, Red7, Red8;
Blue0 = LineP0[-3]; Green0 = LineP0[-2]; Red0 = LineP0[-1];
Blue1 = LineP0[0]; Green1 = LineP0[1]; Red1 = LineP0[2];
Blue2 = LineP0[3]; Green2 = LineP0[4]; Red2 = LineP0[5];
Blue3 = LineP1[-3]; Green3 = LineP1[-2]; Red3 = LineP1[-1];
Blue4 = LineP1[0]; Green4 = LineP1[1]; Red4 = LineP1[2];
Blue5 = LineP1[3]; Green5 = LineP1[4]; Red5 = LineP1[5];
Blue6 = LineP2[-3]; Green6 = LineP2[-2]; Red6 = LineP2[-1];
Blue7 = LineP2[0]; Green7 = LineP2[1]; Red7 = LineP2[2];
Blue8 = LineP2[3]; Green8 = LineP2[4]; Red8 = LineP2[5];
if (Blue1 > Blue2) Swap(Blue1, Blue2);
if (Blue4 > Blue5) Swap(Blue4, Blue5);
if (Blue7 > Blue8) Swap(Blue7, Blue8);
if (Blue0 > Blue1) Swap(Blue0, Blue1);
if (Blue3 > Blue4) Swap(Blue3, Blue4);
if (Blue6 > Blue7) Swap(Blue6, Blue7);
if (Blue1 > Blue2) Swap(Blue1, Blue2);
if (Blue4 > Blue5) Swap(Blue4, Blue5);
if (Blue7 > Blue8) Swap(Blue7, Blue8);
if (Blue0 > Blue3) Swap(Blue0, Blue3);
if (Blue5 > Blue8) Swap(Blue5, Blue8);
if (Blue4 > Blue7) Swap(Blue4, Blue7);
if (Blue3 > Blue6) Swap(Blue3, Blue6);
if (Blue1 > Blue4) Swap(Blue1, Blue4);
if (Blue2 > Blue5) Swap(Blue2, Blue5);
if (Blue4 > Blue7) Swap(Blue4, Blue7);
if (Blue4 > Blue2) Swap(Blue4, Blue2);
if (Blue6 > Blue4) Swap(Blue6, Blue4);
if (Blue4 > Blue2) Swap(Blue4, Blue2);
if (Green1 > Green2) Swap(Green1, Green2);
if (Green4 > Green5) Swap(Green4, Green5);
if (Green7 > Green8) Swap(Green7, Green8);
if (Green0 > Green1) Swap(Green0, Green1);
if (Green3 > Green4) Swap(Green3, Green4);
if (Green6 > Green7) Swap(Green6, Green7);
if (Green1 > Green2) Swap(Green1, Green2);
if (Green4 > Green5) Swap(Green4, Green5);
if (Green7 > Green8) Swap(Green7, Green8);
if (Green0 > Green3) Swap(Green0, Green3);
if (Green5 > Green8) Swap(Green5, Green8);
if (Green4 > Green7) Swap(Green4, Green7);
if (Green3 > Green6) Swap(Green3, Green6);
if (Green1 > Green4) Swap(Green1, Green4);
if (Green2 > Green5) Swap(Green2, Green5);
if (Green4 > Green7) Swap(Green4, Green7);
if (Green4 > Green2) Swap(Green4, Green2);
if (Green6 > Green4) Swap(Green6, Green4);
if (Green4 > Green2) Swap(Green4, Green2);
if (Red1 > Red2) Swap(Red1, Red2);
if (Red4 > Red5) Swap(Red4, Red5);
if (Red7 > Red8) Swap(Red7, Red8);
if (Red0 > Red1) Swap(Red0, Red1);
if (Red3 > Red4) Swap(Red3, Red4);
if (Red6 > Red7) Swap(Red6, Red7);
if (Red1 > Red2) Swap(Red1, Red2);
if (Red4 > Red5) Swap(Red4, Red5);
if (Red7 > Red8) Swap(Red7, Red8);
if (Red0 > Red3) Swap(Red0, Red3);
if (Red5 > Red8) Swap(Red5, Red8);
if (Red4 > Red7) Swap(Red4, Red7);
if (Red3 > Red6) Swap(Red3, Red6);
if (Red1 > Red4) Swap(Red1, Red4);
if (Red2 > Red5) Swap(Red2, Red5);
if (Red4 > Red7) Swap(Red4, Red7);
if (Red4 > Red2) Swap(Red4, Red2);
if (Red6 > Red4) Swap(Red6, Red4);
if (Red4 > Red2) Swap(Red4, Red2);
LinePD[0] = Blue4;
LinePD[1] = Green4;
LinePD[2] = Red4;
LineP0 += 3;
LineP1 += 3;
LineP2 += 3;
LinePD += 3;
}
}
}
}
看上去代码的行数多了，但是实际上执行速度会更快，我们测试的结果如下：
1920*1080大小的24位图像，平均用时155ms，灰度图像平均用时45ms，比之前的原始实现速度要快了近10倍。
而在任意半径中值滤波（扩展至百分比滤波器）O(1)时间复杂度算法的原理、实现及效果一文中的算法，采用了SSE优化，同样大小的图耗时为：
1920*1080大小的24位图像，平均用时260ms，灰度图像平均用时160ms，比上述的C语言版本要慢。
早期有朋友曾提示我在手机上使用Neon可以做到16MB的图像半径为1的中值模糊可以做到20ms，我真的一点也不敢相信。总觉得不太可思议。16MB可是4000*4000的大小啊，我用上述C的代码处理起来要242ms，比手机端还慢了10倍。
经过朋友提醒，在https://github.com/ARM-software/ComputeLibrary/blob/master/src/core/NEON/kernels/NEMedian3x3Kernel.cpp#L113上看到了相关的Neon代码，惊奇的发现他和我上面的C代码几乎完全一样。但是就是这一点代码提醒了我。
inline void sort(uint8x8_t &a, uint8x8_t &b)
{
const uint8x8_t min = vmin_u8(a, b);
const uint8x8_t max = vmax_u8(a, b);
a = min;
b = max;
}
真是一语惊醒梦中人啊，这么简单的优化我怎么没想到呢。
我们自己看看上面的C代码，每个像素的9次比较虽然不能用SIMD指令做，但是多个像素的比较之间是相互不关联的，因此，这样我就可以一次性处理16个像素了，改成SSE优化的方式也就很简单了：
inline void _mm_sort_ab(__m128i &a, __m128i &b)
{
const __m128i min = _mm_min_epu8(a, b);
const __m128i max = _mm_max_epu8(a, b);
a = min; b = max;
}
void MedianBlur3X3_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
int Channel = Stride / Width;
int BlockSize = 16, Block = ((Width - 2)* Channel) / BlockSize;
for (int Y = 1; Y < Height - 1; Y++)
{
unsigned char *LineP0 = Src + (Y - 1) * Stride + Channel;
unsigned char *LineP1 = LineP0 + Stride;
unsigned char *LineP2 = LineP1 + Stride;
unsigned char *LinePD = Dest + Y * Stride + Channel;
for (int X = 0; X < Block * BlockSize; X += BlockSize, LineP0 += BlockSize, LineP1 += BlockSize, LineP2 += BlockSize, LinePD += BlockSize)
{
__m128i P0 = _mm_loadu_si128((__m128i *)(LineP0 - Channel));
__m128i P1 = _mm_loadu_si128((__m128i *)(LineP0 - 0));
__m128i P2 = _mm_loadu_si128((__m128i *)(LineP0 + Channel));
__m128i P3 = _mm_loadu_si128((__m128i *)(LineP1 - Channel));
__m128i P4 = _mm_loadu_si128((__m128i *)(LineP1 - 0));
__m128i P5 = _mm_loadu_si128((__m128i *)(LineP1 + Channel));
__m128i P6 = _mm_loadu_si128((__m128i *)(LineP2 - Channel));
__m128i P7 = _mm_loadu_si128((__m128i *)(LineP2 - 0));
__m128i P8 = _mm_loadu_si128((__m128i *)(LineP2 + Channel));
_mm_sort_ab(P1, P2); _mm_sort_ab(P4, P5); _mm_sort_ab(P7, P8);
_mm_sort_ab(P0, P1); _mm_sort_ab(P3, P4); _mm_sort_ab(P6, P7);
_mm_sort_ab(P1, P2); _mm_sort_ab(P4, P5); _mm_sort_ab(P7, P8);
_mm_sort_ab(P0, P3); _mm_sort_ab(P5, P8); _mm_sort_ab(P4, P7);
_mm_sort_ab(P3, P6); _mm_sort_ab(P1, P4); _mm_sort_ab(P2, P5);
_mm_sort_ab(P4, P7); _mm_sort_ab(P4, P2); _mm_sort_ab(P6, P4);
_mm_sort_ab(P4, P2);
_mm_storeu_si128((__m128i *)LinePD, P4);
}
for (int X = Block * BlockSize; X < (Width - 2) * Channel; X++, LinePD++)
{
// DO Something
}
}
}
注意到上面我已经把灰度和彩色图像的代码写成同一个方式处理了，这是因为对于彩色图像，3个通道之间的处理时毫无联系的。同时我们前面的Swap2个变量的过程时完全可以通过Min和Max两个算子实现的，我们按下F5测试运行，惊人的速度出现了：
1920*1080大小的24位图像，平均用时3ms，灰度图像平均用时1ms，比上述的C语言版本快了近40倍。
顺便也测试了下16MB的图像，结果平均只需要7.5ms。真是太厉害了。
二、5*5中值模糊
对于5*5的中值模糊，优化的方式还是一样的，但是5*5共计25个像素，理论上需要131次比较，其他的过程类似，测试基于SSE的方式，5*5的中值1920*1080大小的24位图像，平均用时40ms，灰度图像平均用时20ms，虽慢了很多，但是还是O(1)那里的速度快。
三、蒙尘和划痕
在这里提Photoshop的这个算法也许并不是很合适，但是我也是在研究中值模糊时顺便把这个算法给攻破的，当我们打开蒙尘和划痕界面时，发现其有半径和阈值两个参数，细心比较，如果阈值设置为0，则相同半径设置时其结果图像和杂色里的中间值算法的结果一模一样，这也可以从蒙尘和划痕算法和中间值同样都放在杂色菜单下可以看出端倪。
通过上述分析，我们可以肯定蒙尘和划痕算法是基于中值模糊的，实际上，PS里很多算法都是基于中值模糊的，特别是那些有平滑度参数的算法^_^。经过多次测试，我们得到的该算法的结果就是如下：
if Abs(Median - Src) > Threshold
Dest = Median
else
Dest = Src
对于彩色图像，不是用彩色图的中值，而是用其亮度值作为唯一的判断标准，如果用彩色的中值作为标准来判断每个分量的，很容易出现过多的噪点，因为有可能会出现Blue分量改变，而Red不变的情况，或其他类似现象。
蒙尘和划痕的一个作用是去除噪点，特别的，我觉得他在小半径的时候更为有用，小半径中值不会改变原图太多，加上这个阈值则可以很容易去除噪点，同时，基本不会出现新的模糊问题。比如下面这个图。
原图 半径为1的中值模糊
半径为1，阈值取20时的蒙尘和划痕 半径为2，阈值取20时的蒙尘和划痕
由以上几图，可以明显的看出，带阈值的蒙尘和划痕在抑制了噪音的同时对原图其他细节基本没有破坏，因此，是一种比较合适的初级的预处理算法，既然是预处理，那么其效率就非常重要了，因此本文的快速3*3模糊的作用也就是相当有用。
还举个例子，下面这个照片中有很多白色的小点点，如果直接用中值确实可以将白点去除，但是可能要半径为四左右才可以去除，但是此时图像整体也变得模糊了，如果使用蒙版和划痕，则处理后的效果非常完美。
原图 半径4，阈值100
本文相关算法代码下载地址：https://files.cnblogs.com/files/Imageshop/MedianBlur3X3.rar。
本人的SSE算法优化合集DEMO:测试Demo:https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar。
