[快速閱讀九] 自適應中值濾波及保守濾波在去除椒鹽噪音或脈衝噪音上的應用。

Imageshop發表於2024-11-04

　　這兩個濾波器也是很久前就看過的，最近偶然翻起那本比較經典的matlab數字影像處理（岡薩雷斯）書，裡面也提到了這個演算法，覺得效果還行，於是也還是稍微整理下。

　　為了自己隨時翻資料舒服和省事省時，這個演算法的原理我們還是把他從別人的部落格裡搬過來吧：

　　摘自：影像處理基礎(2)：自適應中值濾波器(基於OpenCV實現)

　　自適應的中值濾波器也需要一個矩形的視窗Sxy，和常規中值濾波器不同的是這個視窗的大小會在濾波處理的過程中進行改變（增大）。需要注意的是，濾波器的輸出是一個畫素值，該值用來替換點(x,y)處的畫素值，點(x,y)是濾波視窗的中心位置。

原理說明

　　過程A的目的是確定當前視窗內得到中值Zmed是否是噪聲。如果Zmin<Zmed<Zmax，則中值Zmed不是噪聲，這時轉到過程B測試，當前視窗的中心位置的畫素Zxy是否是一個噪聲點。如果Zmin<Zxy<Zmax，則Zxy不是一個噪聲，此時濾波器輸出Zxy；如果不滿足上述條件，則可判定Zxy是噪聲，這是輸出中值Zmed（在A中已經判斷出Zmed不是噪聲）。

　　如果在過程A中，得到則Zmed不符合條件Zmin<Zmed<Zmax，則可判斷得到的中值Zmed是一個噪聲。在這種情況下，需要增大濾波器的視窗尺寸，在一個更大的範圍內尋找一個非噪聲點的中值，直到找到一個非噪聲的中值，跳轉到B；或者，視窗的尺寸達到了最大值，這時返回找到的中值，退出。

　　從上面分析可知，噪聲出現的機率較低，自適應中值濾波器可以較快的得出結果，不需要去增加視窗的尺寸；反之，噪聲的出現的機率較高，則需要增大濾波器的視窗尺寸，這也符合種中值濾波器的特點：噪聲點比較多時，需要更大的濾波器視窗尺寸。

　　摘抄完成..............................................

　　這個過程理解起來也不是很困難，影像處理基礎(2)：自適應中值濾波器(基於OpenCV實現) 這個部落格也給出了參考程式碼，不過我很遺憾的告訴大家，這個部落格的效果雖然可以，但是編碼和很多其他的部落格一樣，是存在問題的。

　　核心在這裡：

for (int j = maxSize / 2; j < im1.rows - maxSize / 2; j++)
    {
        for (int i = maxSize / 2; i < im1.cols * im1.channels() - maxSize / 2; i++)
        {
            im1.at<uchar>(j, i) = adaptiveProcess(im1, j, i, minSize, maxSize);
        }
    }

　　他這裡的就求值始終是對同一份圖，這樣後續的處理時涉及到的領域實際上前面的部分已經是被修改的了，不符合真正的原理的。 至於為什麼最後的結果還比較合適，那是因為這裡的領域相關性不是特別強。

　　我這裡藉助於最大值和最小值濾波以及中值濾波，一個簡單的實現如下所示：

/// <summary>
/// 實現影像的自使用中值模糊。更新時間2015.3.11。
/// 參考：Adaptive Median Filtering Seminar Report By: PENG Lei (ID: 03090345)
/// </summary>
/// <param name="Src">需要處理的源影像的資料結構。</param>
/// <param name="Dest">儲存處理後的影像的資料結構。</param>
/// <param name="Radius">濾波的半徑，有效範圍[1,127]。</param>
/// <remarks> 1: 能處理8位灰度和24位及32點陣圖像。</remarks>
/// <remarks> 2: 半徑大於10以後基本沒區別了。</remarks>
int IM_AdaptiveMedianFilter(unsigned char* Src, unsigned char* Dest, int Width, int Height, int Stride, int MinRadius, int MaxRadius)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                    return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                        return IM_STATUS_INVALIDPARAMETER;
    if ((MinRadius <= 0) || (MaxRadius <= 0))                return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3))                    return IM_STATUS_NOTSUPPORTED;
    int Status = IM_STATUS_OK;
    
    int Threshold = 10;
    if (MinRadius > MaxRadius)    IM_Swap(MinRadius, MaxRadius);

    bool AllProcessed = false;
    unsigned char* MinValue = (unsigned char*)malloc(Height * Stride * sizeof(unsigned char));
    unsigned char* MaxValue = (unsigned char*)malloc(Height * Stride * sizeof(unsigned char));
    unsigned char* MedianValue = (unsigned char*)malloc(Height * Stride * sizeof(unsigned char));
    unsigned char* Flag = (unsigned char*)malloc(Height * Width * sizeof(unsigned char));
    if ((MinValue == NULL) || (MaxValue == NULL) || (MedianValue == NULL) || (Flag == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    memset(Flag, 0, Height * Width * sizeof(unsigned char));

    if (Channel == 1)
    {
        //    The median filter starts at size 3-by-3 and iterates up to size MaxRadius-by-MaxRadius
        for (int Z = MinRadius; Z <= MaxRadius; Z++)
        {
            Status = IM_MinFilter(Src, MinValue, Width, Height, Stride, Z);
            if (Status != IM_STATUS_OK)        goto FreeMemory;
            Status = IM_MaxFilter(Src, MaxValue, Width, Height, Stride, Z);
            if (Status != IM_STATUS_OK)        goto FreeMemory;        
            Status = IM_MedianBlur(Src, MedianValue, Width, Height, Stride, Z, 50);
            if (Status != IM_STATUS_OK)        goto FreeMemory;
        
            for (int Y = 0; Y < Height; Y++)
            {
                int Index = Y * Stride;
                int Pos = Y * Width;
                for (int X = 0; X < Width; X++, Index++, Pos++)
                {
                    if (Flag[Pos] == 0)
                    {
                        int Min = MinValue[Index], Max = MaxValue[Index], Median = MedianValue[Index];
                        if ((Median > Min + Threshold) && (Median < Max - Threshold))
                        {
                            int Value = Src[Index];
                            if ((Value > Min + Threshold) && (Value < Max- Threshold))
                            {
                                Dest[Index] = Value;
                            }
                            else
                            {
                                Dest[Index] = Median;
                            }
                            Flag[Pos] = 1;
                        }
                    }
                }
            }
            AllProcessed = true;
            for (int Y = 0; Y < Height; Y++)
            {
                int Pos = Y * Width;
                for (int X = 0; X < Width; X++)
                {
                    if (Flag[Pos + X] == 0)        //    檢測是否每個點都已經處理好了
                    {
                        AllProcessed = false;
                        break;
                    }
                }
                if (AllProcessed == false) break;
            }
            if (AllProcessed == true) break;
        }

        /*    Output zmed for any remaining unprocessed pixels. Note that this
            zmed was computed using a window of size Smax-by-Smax, which is
            the final value of k in the loop.*/

        if (AllProcessed == false)
        {
            for (int Y = 0; Y < Height; Y++)
            {
                int Index = Y * Stride;
                int Pos = Y * Width;
                for (int X = 0; X < Width; X++)
                {
                    if (Flag[Pos + X] == 0)    Dest[Index + X] = Src[Index + X];
                }
            }
        }
    }
    else
    {

    }
FreeMemory:
    if (MinValue != NULL)        free(MinValue);
    if (MaxValue != NULL)        free(MaxValue);
    if (MedianValue != NULL)    free(MedianValue);
    if (Flag != NULL)            free(Flag);
    return Status;
}

　　注意這裡，我們還做了適當的修改，增加了一個控制閾值Threshold，把原先的　if ((Median > Min) && (Median < Max))

修改為：

　　if ((Median > Min + Threshold) && (Median < Max - Threshold))

　　也可以說是對應用場景的一種擴充套件，增加了函式的韌性。

　　當我們的噪音比較少時，這個函式會很快收斂，也就是不需要進行多次的計算的。

　　另外還有一個演算法叫Conservative Smoothing，翻譯成中文可以稱之為保守濾波，這個演算法在https://homepages.inf.ed.ac.uk/rbf/HIPR2/csmooth.htm有較為詳細的說明，其基本原理是：

　This is accomplished by a procedure which first finds the minimum and maximum intensity values of all the pixels within a windowed region around the pixel in question. If the intensity of the central pixel lies within the intensity range
spread of its neighbors, it is passed on to the output image unchanged. However, if the central pixel intensity is greater than the maximum value, it is set equal to the maximum value; if the central pixel intensity is less than the minimum value,
it is set equal to the minimum value. Figure 1 illustrates this idea.

　　

　　比如上圖的3*3領域，除去中心點之外的8個點其最大值為127，最小值是115，而中心點的值150大於這個最大值，所以中心點的值會被修改為127。

　　注意這裡和前面的自適應中值濾波有一些不同。

　　1、雖然他也利用到了領域的最大值和最小值，但是這個領域是不包含中心畫素本身的，這個和自適應中值是最大的區別。

2、這個演算法在滿足某個條件時，不是用中值代替原始畫素，而是用前面的改造的最大值或者最小值。

　　3、這個演算法也可以改造成和自適應中值一樣，不斷的擴大半徑。

　　4、對於同一個半徑，這個函式多次迭代效果不會有區別。

　　一個簡單的實現如下：

　　　　for (int X = 0; X < Width * Channel; X++, LinePD++)
        {
            int P0 = First[X], P1 = First[X + Channel], P2 = First[X + 2 * Channel];
            int P3 = Second[X], P4 = Second[X + Channel], P5 = Second[X + 2 * Channel];
            int P6 = Third[X], P7 = Third[X + Channel], P8 = Third[X + 2 * Channel];

            int Max0 = IM_Max(P0, P1);
            int Min0 = IM_Min(P0, P1);

            int Max1 = IM_Max(P2, P3);
            int Min1 = IM_Min(P2, P3);

            int Max2 = IM_Max(P5, P6);
            int Min2 = IM_Min(P5, P6);

            int Max3 = IM_Max(P7, P8);
            int Min3 = IM_Min(P7, P8);

            int Max = IM_Max(IM_Max(Max0, Max1), IM_Max(Max2, Max3));
            int Min = IM_Min(IM_Min(Min0, Min1), IM_Min(Min2, Min3));

            if (P4 > Max)
                P4 = Max;
            else if (P4 < Min)
                P4 = Min;
            LinePD[0] = P4;
        }

　　因為去除了中心點後進行的最大值和最小計算，這個演算法如果要實現高效率的和半徑無關的版本，還是需要做一番研究的，目前我尚沒有做出成果。

　　我們拿標準的Lena圖測試，使用matlab分別為他增加0.02及0.2範圍的椒鹽噪音，然後使用自適應中值濾波及保守濾波進行處理。

　　　　　　　　　　新增0.02機率的椒鹽噪音圖　　　　　　　　　　　　　　　　　　　　　　　 　最小半徑1，最大半徑13時的去噪效果　　　　　　　　　　　　　　　　　　　　　　　　　　3*3的保守濾波

　　　　　　新增0.2機率的椒鹽噪音圖　　　　　　　　　　　　　　　　　　　　　　　　　　　 　　最小半徑1，最大半徑5時的去噪效果　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　3*3的保守濾波

　　可以看到，自適應中值濾波在去除椒鹽噪音的效果簡直就是逆天，基本完美的復現了原圖，有的時候我自己都不敢相信這個結果。而保守濾波由於目前我只實現了3*3的版本，因此對於噪音比較集中的地方暫時無法去除，這個可以透過擴大半徑和類似使用自適應中值的方式處理，而這種噪音的影像使用DCT去噪、小波去噪、NLM等等都不具有很好的效果，唯一能和他比擬的就只有蒙版和劃痕裡面小引數時的效果（這也是以中值為基礎的）。所以去除椒鹽噪音還是得靠中值相關的演算法啊。

　　和常規的中值濾波器相比，自適應中值濾波器能夠更好的保護影像中的邊緣細節部分，當然代價就是增加了演算法的時間複雜度。

相關文章