<<一種基於δ函數的圖象邊緣檢測算法>>一文算法的實現。

　　原始論文下載： 一種基於δ函數的圖象邊緣檢測算法。

     這篇論文讀起來感覺不像現在的很多論文，廢話一大堆，而是直入主題，反倒使人覺得文章的前后跳躍有點大，不過算法的原理已經講的清晰了。

    一、原理

     文中提出的邊緣檢測算法原理也不是特別復雜，使用了一個低通濾波函數以及一個高通濾波函數，其形式分別為：

                                                     （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#的程序：http://files.cnblogs.com/Imageshop/EdgeDetectUseDeltaFunction.rar

 

*********************************作者： laviewpbt   時間： 2013.10.26    聯系QQ:  33184777  轉載請保留本行信息************************