前言:
而今同学们对“均值滤波算法公式”大体比较关怀,各位老铁们都想要知道一些“均值滤波算法公式”的相关文章。那么小编也在网上网罗了一些对于“均值滤波算法公式””的相关资讯,希望小伙伴们能喜欢,咱们快快来学习一下吧!关于高斯模糊,我在我早期的博客里也有两篇文章予以描述:
SSE图像算法优化系列二:高斯模糊算法的全面优化过程分享(一)。
SSE图像算法优化系列二:高斯模糊算法的全面优化过程分享(二)。
一个是递归的IIR滤波器,一个Deriche滤波器,他们的速度都已经是顶级的了,而且都能够使用SIMD指令优化,其中有讲到《Recursive implementation of the Gaussian filter》这个方法在半径较大的时候会出现一定的瑕疵,核心原因是大半径会导致其中的某些系数特别小,因此造成浮点精度的丢失,因此,要保证效果就必须在计算过程中使用double数据类型,而使用了double,普通的sse指令集的增速效果就不是很明显了,因此,为了速度可能需要使用AVX或者更高的AVX512。
那么这两个都是从离散角度来说比较精确的算法,因为有了SIMD指令,使得他们在PC上即使有了大量的浮点计算,计算的速度也比较不错。在一些特殊的情况或者特殊的硬件中,还是存在浮点计算非常慢的情况(比如FPGA),因此,还是有不上资料和文章提出了一些近似的算法来减少计算量,在celerychen分享的文章快速高斯滤镜算法中,除了上述两个算法外,还提到了均值滤波逼近高斯滤波以及 扩展二项式滤波逼近高斯滤波两个方法。
一、Binomial Filter 二项式滤波滤波器
多年前我也看过这个文章,那个时候也没有怎么在意,最近在研究halcon的一些滤波器时,偶尔翻到其binomial_filter函数的说明时,突然看到其有如下的一些描述:
The binomial filter is a very good approximation of a Gaussian filter that can be implemented extremely efficiently using only integer operations. Hence, binomial_filter is very fast.
他说这个binomial filter只有整形的计算,因此速度非常快。
是不是这样的呢,于是我在网络上搜索这方面的资料,大部分都指向 extended binomial filter for Fast gaussian Blur这篇文章,可以参考道客巴巴里的这篇文章:
这个文章写得很详细,有着非常丰富的数学推导和公式,还有很多漂亮的曲线,比如下面这个描述了不同的半径r和不同的n阶毕竟的精度。当n=3或者n=4时,也已经相当的准确了。特别是半径比价大时(模糊程度)。
这个文章的末尾提供了相关的参考代码,但是那个参考代码有点凌乱,而且不知道是啥语言,好像是可以为PS写插件的脚本一样,不过大概还是能猜到是什么意思。
celerychen在他的文章说里面几个重要的初始参数不知道如何确定,主要可能是指的scaleFactor,这个其实因该他做界面预览时的一些缩放值,我们做算法时,他就是1。
另外,那个文章是比较早的了,作者提供的那个网站里有提供最新的一批代码,详细见:Efficient Gaussian Blur Algorithm
这个文章中提供的一个重要的计算公式即为:
其中
为用户输入的标准差,n为离散的阶数,r为计算得到半径,由上式我们可以推到出:
r = sqrtf(12.0 * Sigma * Sigma / n + 1)
注意,这里得到的r应该是浮点数,但是浮点数,我们无法进行均值模糊,所以一般需要四舍五入取整。当然,如果要求精度的,那就要去上下两个半径值分别做处理后,在对结果进行插值。
这个公式在 均值滤波逼近高斯滤波 的文章里也有提到。
我们从Efficient Gaussian Blur Algorithm 这个网站的相关链接里可以有如下描述:
gauss2.ffp: The program for the first degree approximation is a great improvement to the box blur. It is as easy to understand and has only slight approximation error.gauss3.ffp: The program for the second degree approximation is accurate enough for most applications.gauss4.ffp: The third degree approximation has nearly no difference to the Gaussian blur for 8-bit images apart of rounding errors
我想说的是gauss2.ffp其实是Third degree approximation对应的代码(n = 3),gauss3.ffp和gauss4.ffp都是Fourth degree approximation对应的代码(n=4),但是他们里面的叠加顺序写法又有所不同(里面的变量sum, deri1, deri2, diff叠加的顺序),那么真正的Second degree approximation文章里没有给链接,但是可以用在链接 里打开。感觉像是作者故意隐藏了一样。
另外还要注意在论文里的sum, deri1, deri2, diff和链接里的顺序也不太一样。
为了方便大家测试,我从链接里整理了一个可以运行的代码,很丑陋,但是确实效果和速度都还不错。
int IM_BinomialFilter4_PureC_Gray_Raw(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Sigma){ int Channel = Stride / Width; if ((Src == NULL) || (Dest == NULL) || (Dest == NULL)) return IM_STATUS_NULLREFRENCE; if ((Width <= 0) || (Height <= 0)) return IM_STATUS_INVALIDPARAMETER; if (Channel != 1) return IM_STATUS_NOTSUPPORTED; int Status = IM_STATUS_OK; Sigma = IM_Min(Sigma, 256); // 如果全部用整形的数据,半径不能大于30, Sum和Deri用浮点半径可以取得较大 int Radius = (int)(sqrtf(3.0 * Sigma * Sigma + 1) + 0.4999999f); int Weight = Radius * Radius * Radius * Radius; float InvWeight = 1.0f / Weight; unsigned char *Buffer = (unsigned char *)malloc((IM_Max(Width, Height) + 4 * Radius) * sizeof(unsigned char)); unsigned char *Temp = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char)); int *RowOffset = (int *)malloc((Width + 4 * Radius) * sizeof(int)); int *ColOffset = (int *)malloc((Height + 4 * Radius) * sizeof(int)); if ((Buffer == NULL) || (Temp == NULL) || (RowOffset == NULL) || (ColOffset == NULL)) { Status = IM_STATUS_OUTOFMEMORY; goto FreeMemory; } // 注意起点位置是偏移了一个1个像素的 for (int X = 0; X < Width + 4 * Radius; X++) RowOffset[X] = IM_GetMirrorPos(Width, X - 2 * Radius + 1); for (int Y = 0; Y < Height + 4 * Radius; Y++) ColOffset[Y] = IM_GetMirrorPos(Height, Y - 2 * Radius + 1); for (int Y = 0; Y < Height; Y++) // horizontal blur... { unsigned char *LinePS = Src + Y * Stride; unsigned char *LinePD = Temp + Y * Stride; for (int X = 0; X < 2 * Radius - 1; X++) Buffer[X] = LinePS[RowOffset[X]]; memcpy(Buffer + 2 * Radius - 1, LinePS, Width); for (int X = 2 * Radius - 1 + Width; X < Width + 4 * Radius; X++) Buffer[X] = LinePS[RowOffset[X]]; float Sum = 0, Deri1 = 0; int Deri2 = 0, Diff = 0; for (int X = -4 * Radius; X < 0; X++) { Sum += Deri1; Deri1 += Deri2; Deri2 += Diff; // accumulate pixel blur Diff += Buffer[X + 4 * Radius]; if (X + 3 * Radius >= 0) Diff -= 4 * Buffer[X + 3 * Radius]; // -4, if (X + 2 * Radius >= 0) Diff += 6 * Buffer[X + 2 * Radius]; // +6, if (X + Radius >= 0) Diff -= 4 * Buffer[X + 1 * Radius]; // -4, (+1)} } for (int X = 0; X < Width; X++) { Sum += Deri1; Deri1 += Deri2; Deri2 += Diff; Diff += Buffer[X] - 4 * Buffer[X + Radius] + 6 * Buffer[X + 2 * Radius] - 4 * Buffer[X + 3 * Radius] + Buffer[X + 4 * Radius]; //{+1, -4, +6, -4, +1} LinePD[X] = (int)(Sum * InvWeight + 0.499999f); //LinePD[X] = (Sum + Weight / 2) / Weight; } } for (int X = 0; X < Width; X++) // Vertical Blur { unsigned char *LinePS = Temp + X; unsigned char *LinePD = Dest + X; for (int Y = 0; Y < 2 * Radius - 1; Y++) Buffer[Y] = LinePS[ColOffset[Y] * Stride]; for (int Y = 0; Y < Height; Y++) { Buffer[Y + 2 * Radius - 1] = LinePS[Y * Stride]; } for (int Y = 2 * Radius - 1 + Height; Y < Height + 4 * Radius; Y++) Buffer[Y] = LinePS[ColOffset[Y] * Stride]; float Sum = 0, Deri1 = 0; int Deri2 = 0, Diff = 0; for (int Y = -4 * Radius; Y < 0; Y++) { Sum += Deri1; Deri1 += Deri2; Deri2 += Diff; // accumulate pixel blur Diff += Buffer[Y + 4 * Radius]; if (Y + 3 * Radius >= 0) Diff -= 4 * Buffer[Y + 3 * Radius]; // -4, if (Y + 2 * Radius >= 0) Diff += 6 * Buffer[Y + 2 * Radius]; // +6, if (Y + Radius >= 0) Diff -= 4 * Buffer[Y + 1 * Radius]; // -4, (+1)} } for (int Y = 0; Y < Height; Y++) { Sum += Deri1; Deri1 += Deri2; Deri2 += Diff; Diff += Buffer[Y] - 4 * Buffer[Y + Radius] + 6 * Buffer[Y + 2 * Radius] - 4 * Buffer[Y + 3 * Radius] + Buffer[Y + 4 * Radius]; //{+1, -4, +6, -4, +1} LinePD[0] = (int)(Sum * InvWeight + 0.499999f); //LinePD[0] = (Sum + Weight / 2) / Weight; LinePD += Stride; } }FreeMemory: if (Buffer != NULL) free(Buffer); if (Temp != NULL) free(Temp); if (RowOffset != NULL) free(RowOffset); if (ColOffset != NULL) free(ColOffset); return Status;}
注意这个是翻译自gauss3.ffp的代码并做了适当的调整,我们看到程序还是使用了很多的浮点数的,否则当Sigma稍微大一点,结果就溢出了,当然,代码里也给出了如果SIgma不大于30,则可以全部用整形实现。
我们处理了一副3000*2000的灰度图,出乎意料的是,就这个很随意的代码能获得82毫秒的成绩。
优化:
乍一看起来,Vertical Blur那一部分的代码访问内存的方式很不友好,有着大量的跨行访问内存的地方,确实原始的代码是有这个问题,但是这种模糊都有这个特性,在 高斯模糊算法的全面优化过程分享(一) 中我们探讨过垂直方向处理算法一般不宜直接写,而应该用一个临时的行缓存进行处理,这里同样可以用这个方法来处理,而且处理的代码还能通用各种不同的通道,这里不多赘述。
Horizontal blur部分,看代码知道是存在着严重的像素前后依赖的,这种特性非常不利于指令集优化,处理的方式呢,我在博客里也多次提到,把多行(比如四行)合并到一行后,进行处理,这样合并的一行里就有连续4个像素的处理是毫无关系的,就可以借助于SSE之类的指令予以处理了,这部分方式可以参考我的博文: SSE图像算法优化系列五:超高速指数模糊算法的实现和优化(10000*10000在100ms左右实现)。
经过一番折腾呢,这个算法可以优化到13ms左右(4阶),当然,这个代码明显有个特性,随着Sigma的增大,程序里的边缘部分的计算量也会增加,因此,算法的耗时和参数有一定的关系,但是不是特别明显。
如果是使用的2阶算法,则内部计算要更加的少,同样的图片大约在11ms左右可以处理完成,而原始的Deriche算法嘛,耗时大概在18ms上下。
2阶的有个比较好的特性,就是他的全部的计算过程里的累加量都在整形的有效范围呢,因此确实可以用整形计算实现效果。
完成了这个算法后,回头在去看Halcon的binomial_filter ,感觉味道又变了,本来以为他的那个MaskWidth和MaskHeight就是对应的水平和垂直方向的Sigma的,结果测试不是,他的结果在MaskWidth和MaskHeight很大时,模糊的程度也比较小。目前还搞不清楚是怎么回事。
不过有一点,Halcon的这个算子的速度那真叫一个牛逼,3000*2000的图都只要几豪秒(随MaskWidth和MaskHegiht的增加耗时也在增加)。
二、使用多次均值模糊模拟高斯模糊
这个算法的参考文献的正式名字应该是Fast Almost-Gaussian Filtering,而不是celerychen文章里提的 Arbitrary Gaussian Filtering with 25 Addtions and 5 Multiplications per Pixel,那个文章里根本没看到说25次加法和5次乘法。
这个文章的理论基础其实也是这个公式:
英文的原文为:
他这里给出了用均值模糊模拟高斯模糊的半径的计算方式,即上图中公式5,当给定均方差,给定均值模拟的次数,就可以计算出m。
在这里作者给出了相关的简易的示例代码,稍作整理如下所示:
// 使用 N 个 BoxBlur 拟合 Sigma的高斯函数// 参考:Arbitrary Gaussian Filtering with 25 Addtions and 5 Multiplications per Pixel.pdf// Fast Almost - Gaussian Filtering.pdfvoid IM_GetBoxRadiusForGauss(float Sigma, int *Size, int N){ float IdealW = sqrt(12.0 * Sigma * Sigma / N + 1.0); // 论文公式1 int LowerW = floor(IdealW); // 向下取整 if (LowerW % 2 == 0) LowerW--; // the first odd integer that is less than IdealW int UpperW = LowerW + 2; // the other being being the first odd integer greater than IdealW int M = round((12.0 * Sigma * Sigma - N * LowerW * LowerW - 4 * N * LowerW - 3 * N) / (-4 * LowerW - 4)); // 论文公式2 for (int Y = 0; Y < N; Y++) // Apply an averaging filter of width LowerW M times. Size[Y] = (Y < M ? LowerW : UpperW); // Apply an averaging filter of width UpperW N - M times.}int IM_FastGaussBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius){ int Channel = Stride / Width; if ((Src == NULL) || (Dest == NULL) || (Dest == NULL)) return IM_STATUS_NULLREFRENCE; if ((Width <= 0) || (Height <= 0)) return IM_STATUS_INVALIDPARAMETER; if ((Channel != 1) && (Channel != 3) && (Channel != 4)) return IM_STATUS_NOTSUPPORTED; int Status = IM_STATUS_OK; unsigned char *Temp = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char)); if (Temp == NULL) { Status = IM_STATUS_OUTOFMEMORY; goto FreeMemory; } int BoxRadius[3]; // 使用3次均值模糊代替高斯 IM_GetBoxRadiusForGauss(Radius, BoxRadius, 3); // 计算拟合的半径 Status = IM_BoxBlur(Src, Dest, Width, Height, Stride, (BoxRadius[0] - 1) / 2); if (Status != IM_STATUS_OK) goto FreeMemory; Status = IM_BoxBlur(Dest, Temp, Width, Height, Stride, (BoxRadius[1] - 1) / 2); if (Status != IM_STATUS_OK) goto FreeMemory; Status = IM_BoxBlur(Temp, Dest, Width, Height, Stride, (BoxRadius[2] - 1) / 2); if (Status != IM_STATUS_OK) goto FreeMemory;FreeMemory: if (Temp != NULL) free(Temp); return Status;}
可以看到,就是简单的3次均值模糊,因为均值模糊可以高度的向量化编程,拥有特别特别快的速度,因此,也能获得非常优异的性能。
三、效果比较
对标准的高斯模糊,二阶二项式、4阶二项式以及均值模糊模拟进行测试,发现他们在视觉上无特备明显的差异。
原始图像 Deriche滤波器 2阶多项式滤波器
4阶多项式滤波器 三级均值模糊
可以通过我的DEMO里动态的看到这些滤波器的结果差异,感觉这些差异是在可接收的范围内的。
文章来自
标签: #均值滤波算法公式