<<一种基于δ函数的图象边缘检测算法>>一文算法的实现。
原始论文下载: 一种基于δ函数的图象边缘检测算法。
这篇论文读起来感觉不像现在的很多论文,废话一大堆,而是直入主题,反倒使人觉得文章的前后跳跃有点大,不过算法的原理已经讲的清晰了。
一、原理
文中提出的边缘检测算法原理也不是特别复杂,使用了一个低通滤波函数以及一个高通滤波函数,其形式分别为:
(1)
(2)
当图像中的噪音比较少时,可以直接使用高通滤波器对图像进行滤波,得到图像的细节信息(即边缘处),论文中称之为D算法,计算公式如下:
式中顶部的横线应该是表示开平方的意思。
而当图像含有噪音时,则采用高通和低通滤波器结合方式,使用低通滤波器平滑图像中的噪音,高通滤波器检测边缘,这个原理则类似于高斯拉普拉斯边缘检测过程,论文中称之为C算法,计算公式如下:
式中w表示的是窗口大小,取值越大,边缘的宽度越大,建议理想取值为2。
上面两个式子都已经是离散化的表达方式了,因此实际上也是一种对图像的模板操作,只是模板中的因子需要随着参数的不同而改变。
注意:D算法仅仅是一维的模板操作,而C算法是二维的。
二、代码
下面贴出D算法的核心代码:
void EdgeDetail(byte* Src, byte* Dest, int Width, int Height, int Stride, int Radius = 2, double S = 1, double T = 3) { int X, Y, I, J, XX, YY; byte* SrcP, DestP; int SumOne, SumTwo, Power; byte* SqrValue = (byte*)GlobalAlloc(GPTR, (256 * 256) * sizeof(byte)); int* SpeedHigh = (int*)GlobalAlloc(GPTR, (Radius * 2 + 1) * sizeof(int)); SpeedHigh += Radius; for (Y = 0; Y < 256 * 256; Y++) SqrValue[Y] = (byte)Math.Sqrt(Y); for (Y = -Radius; Y <= Radius; Y++) { if (Y == 0) SpeedHigh[Y] = 0; else SpeedHigh[Y] = (int)((((Math.Cos(S * Y) / Y) - (Math.Sin(S * Y) / S) * (1.0 / (Y * Y) + 1.0 / (T * T))) * Math.Exp(-((double)Y * Y) / (2 * T * T))) * 1024); } for (Y = 0; Y < Height; Y++) { DestP = Dest + Y * Stride; for (X = 0; X < Width; X++) { SumOne = 0; SumTwo = 0; for (J = -Radius; J <= Radius; J++) { XX = X + J; if (XX < 0) XX = 0; else if (XX >= Width) XX = Width - 1; SrcP = Src + Stride * Y + XX; SumOne += (SpeedHigh[J] * SrcP[0]) >> 10; YY = Y + J; if (YY < 0) YY = 0; else if (YY >= Height) YY = Height - 1; SrcP = Src + Stride * YY + X; SumTwo += (SpeedHigh[J] * SrcP[0]) >> 10; } Power = SumOne * SumOne + SumTwo * SumTwo; if (Power > 65025) Power = 65025; DestP[0] = SqrValue[Power]; DestP++; } } SpeedHigh -= Radius; GlobalFree((IntPtr)SqrValue); GlobalFree((IntPtr)SpeedHigh); }
如上所示,我采用了整数运算代替了浮点运算,主要目的是为了提高速度,当然这样做可能会牺牲一部分精度,由于从算法的必要性上讲,Radius不需要取得很大,因此,对于内部的二重循环来说,压力不是特大,因此没有做特殊的优化。而在超出边界处,直接采用的是使用边界元素值。
上述代码的内部循环里有一些计算式可以提取到外部来的, 只是为了算法的清晰性,未做优化,速度发烧友可以自行提取。
该算法各像素之间的计算式独立的,因此可以很简单的就实现并行计算。
而C算法的代码就稍微复杂一点:
void EdgeCoarse(byte* Src, byte* Dest, int Width, int Height, int Stride, int Radius = 2, double S0 = 0.3, double T0 = 3, double S1 = 0.2, double T1 = 2) { int X, Y, I, J, XX, YY; byte* SrcP, DestP; int SumOne, SumTwo, Power; int* SqrValue = (int*)GlobalAlloc(GPTR, (256 * 256) * sizeof(int)); int* SpeedHigh = (int*)GlobalAlloc(GPTR, (Radius * 2 + 1) * sizeof(int)); int* SpeedLow = (int*)GlobalAlloc(GPTR, (Radius * 2 + 1) * sizeof(int)); SpeedHigh += Radius; SpeedLow += Radius; for (Y = 0; Y < 256 * 256; Y++) SqrValue[Y] = (int)Math.Sqrt(Y); for (Y = -Radius; Y <= Radius; Y++) { if (Y == 0) { SpeedHigh[Y] = 0; SpeedLow[Y] = 1024; } else { SpeedHigh[Y] = (int)((((Math.Cos(S1 * Y) / Y) - (Math.Sin(S1 * Y) / S1) * (1.0 / (Y * Y) + 1.0 / (T1 * T1))) * Math.Exp(-((double)Y * Y) / (2 * T1 * T1))) * 1024); SpeedLow[Y] = (int)(((Math.Sin(S0 * Y) / (S0 * Y)) * Math.Exp(-((double)Y * Y) / (2 * T0 * T0))) * 1024); } } for (Y = 0; Y < Height; Y++) { DestP = Dest + Y * Stride; for (X = 0; X < Width; X++) { SumOne = 0; SumTwo = 0; for (J = -Radius; J <= Radius; J++) { YY = Y + J; if (YY < 0) YY = 0; else if (YY >= Height) YY = Height - 1; for (I = -Radius; I <= Radius; I++) { XX = X + I; if (XX < 0) XX = 0; else if (XX >= Width) XX = Width - 1; SrcP = Src + Stride * YY + XX; SumOne += (SpeedHigh[I] * SpeedLow[J] * SrcP[0]) >>20; SumTwo += (SpeedLow[I] * SpeedHigh[J] * SrcP[0]) >>20; } } Power = SumOne * SumOne + SumTwo * SumTwo; if (Power > 65025) Power = 65025; DestP[0] = (byte)SqrValue[Power]; DestP++; } } SpeedHigh -= Radius; SpeedLow -= Radius; GlobalFree((IntPtr)SqrValue); GlobalFree((IntPtr)SpeedHigh); GlobalFree((IntPtr)SpeedLow); }
我个人不怎么喜欢用C#的数组,这也是从性能角度考虑的,我喜欢直接操作指针。这个可以根据每个人自己的习惯修改吧。
相信能看懂原理的朋友对于代码部分的理解也应该很容易,这里不做多解释。
三、效果
c算法的结果
原图 Radius=2,S=3.14,T=1 Radius=2,S=1.57,T=1
D算法:
原图 Radius=2,S0 = 0.3, T0 = 3, S1 = 0.2, T1 = 2 Radius=2,S0 = 3, T0 = 3, S1 = 2, T1 = 2
可见,这个算法要取得比较好的效果,是需要调整S/T这些参数,关于这些参数的取值意向,可以参考原文中的一些描述。
这个工程比较简单,附上C#的程序:https://files.cnblogs.com/Imageshop/EdgeDetectUseDeltaFunction.rar
*********************************作者: laviewpbt 时间: 2013.10.26 联系QQ: 33184777 转载请保留本行信息************************