基於value-and-criterion structure方式的實現的濾波器在原理上其實比較簡單，感覺下面論文中得一段話已經描述的比較清晰了，直接貼英文吧，感覺翻譯過來反而失去了原始的韻味了。

The value-and-criterion filter structure is based on the geometrical structure of mathematical morphology, but allows the use of a much wider variety of operations (both linear and nonlinear) than standard morphology. Morphological filters usually only use extreme order statistic operations (maximum and minimum). The new filter structure is much more flexible, but still retains the geometric characteristics of mathematical morphology.

The value-and-criterion filter structure is derived from the “compound” structure of the morphological operations opening and closing. Opening and closing are both sequential operations; opening consists of the successive application of the basic morphological operations of erosion and dilation, and closing is dilation followed by erosion. Like opening or closing, a value- and-criterion function has two distinct stages. However, a value-and-criterion filter can have two separate operations acting in parallel in the first stage, as opposed to a single operation (such as erosion in the case of opening). The second stage of a value-and-criterion function uses the output of one of these first stage operations to determine which output from the other first stage operation is selected as the final output。

　　簡單的理解，value-and-criterion structure這種結構的濾波器有兩個元素，一個是值(Value)函數，一個是評判標準(criterion)函數，普通的過程我們我們只需要計算出值函數，就結束了，而value-and-criterion structure這種結構需要根據評判標準來選擇最後的值函數的值。

　　經典的此類函數有Kuwahara濾波器，Mean of Least Variance(MLV)濾波器，Minimum Coefficient of Variation(MCV)濾波器等，我們稍微介紹一下。

　　一、Kuwahara濾波器

很明顯，這個是由Kuwahara等人提出來的。其基本原理是：計算圖像模板中鄰域內的均值和方差，選擇圖像灰度值較為均勻的區域的均值替代模板中心像素灰度值。而圖像模板中較為均勻的區域所對應的方差是最小的。為了獲取圖像模板中較為均勻區域的均值，濾波區域R被劃分為K個重疊的子區域R1，R2,…, Rk。位於圖像中位置為（u，v）處的每一個像素，其所對應的每一個模板子區域的均值和方差都將被計算。那麽標準的Kuwahara濾波器只計算四塊鄰域的方差，取方差最小的鄰域計算其平均值，得到的結果作為目標像素的新值。當然後來還有Tomita-Tsuji，Nagao-Matsuyama’s等人提出了不同的領域，詳見https://blog.csdn.net/lz0499/article/details/54646952一文。

　　在這個濾波器裏，value函數對應的是均值，criterion 函數對應的是均方差，從算法實現上核心的是均值和方差的快速實現，這在我的博文SSE圖像算法優化系列十四：局部均方差及局部平方差算法的優化 裏已經有描述，此處不贅述。

那麽一個簡單的Kuwahara濾波器的主要代碼可能如下：

int IM_KuwaharaFilter(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                                return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                                    return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3))                                return IM_STATUS_INVALIDPARAMETER; 
    
    int Status = IM_STATUS_OK;
    unsigned char *Average = (unsigned char *)malloc(Height * Width * sizeof(unsigned char));
    int *Deviation = (int *)malloc(Height * Width * sizeof(int));
    int *RowOffset = (int *)malloc((Width + Radius + Radius) * sizeof(int));
    int *ColOffset = (int *)malloc((Height + Radius + Radius) * sizeof(int));
    unsigned char *Blur = NULL;

    if ((Average == NULL) || (Deviation == NULL) || (RowOffset == NULL) || (ColOffset == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    if (Channel == 1)
    {
        Status = IM_GetLocalMeanAndDeviationFilter(Src, Average, Deviation, Width, Height, Stride, Radius, false);
        if (Status != IM_STATUS_OK)    goto FreeMemory;
    }
    else
    {
        Blur = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));        //    彩色需要使用明度通道來計算均方差，不然會出現彩色異常塊
        if (Blur == NULL)
        {
            Status = IM_STATUS_OUTOFMEMORY;
            goto FreeMemory;
        }

        Status = IM_GetLuminance(Src, Average, Width, Height, Stride, false);            //    得到明度通道
        if (Status != IM_STATUS_OK)    goto FreeMemory;
        Status = IM_GetLocalMeanAndDeviationFilter(Average, NULL, Deviation, Width, Height, Width, Radius, false);    //    無需保存對應的均值，因為沒用
        if (Status != IM_STATUS_OK)    goto FreeMemory;
        Status = IM_BoxBlur(Src, Blur, Width, Height, Stride, Radius);                    //    彩色模糊
        if (Status != IM_STATUS_OK)    goto FreeMemory;
    }

    for (int Y = 0; Y < Height + Radius + Radius; Y++)
        ColOffset[Y] = IM_GetMirrorPos(Height, Y - Radius);

    for (int X = 0; X < Width + Radius + Radius; X++)
        RowOffset[X] = IM_GetMirrorPos(Width, X - Radius);

    for (int Y = 0; Y < Height; Y++)
    {
        int *LinePST = Deviation + ColOffset[Y] * Width;            //    均方差
        int *LinePSB = Deviation + ColOffset[Y + Radius + Radius] * Width;
        unsigned char *LinePD = Dest + Y * Stride;
        if (Channel == 1)
        {
            unsigned char *LinePT = Average + ColOffset[Y] * Width;        //    均值
            unsigned char *LinePB = Average + ColOffset[Y + Radius + Radius] * Width;
            for (int X = 0; X < Width; X++)
            {
                int OffsetL = RowOffset[X];                                            //    求4個點的均方差的最小值
                int OffsetR = RowOffset[X + Radius + Radius];
                int Min, Mean;
                if (LinePST[OffsetL] < LinePST[OffsetR])
                {
                    Min = LinePST[OffsetL];
                    Mean = LinePT[OffsetL];
                }
                else
                {
                    Min = LinePST[OffsetR];
                    Mean = LinePT[OffsetR];
                }
                if (Min > LinePSB[OffsetL])
                {
                    Min = LinePSB[OffsetL];
                    Mean = LinePB[OffsetL];
                }
                if (Min > LinePSB[OffsetR])
                {
                    Min = LinePSB[OffsetR];
                    Mean = LinePB[OffsetR];
                }
                LinePD[X] = Mean;
            }
        }
        else
        {
            unsigned char *LinePT = Blur + ColOffset[Y] * Stride;        //    均值
            unsigned char *LinePB = Blur + ColOffset[Y + Radius + Radius] * Stride;
            for (int X = 0; X < Width; X++)
            {
                int OffsetL = RowOffset[X];                                            //    求4個點的均方差的最小值
                int OffsetR = RowOffset[X + Radius + Radius];
                int Min, MeanB, MeanG, MeanR;
                if (LinePST[OffsetL] < LinePST[OffsetR])
                {
                    Min = LinePST[OffsetL];
                    MeanB = LinePT[OffsetL * 3 + 0];
                    MeanG = LinePT[OffsetL * 3 + 1];
                    MeanR = LinePT[OffsetL * 3 + 2];
                }
                else
                {
                    Min = LinePST[OffsetR];
                    MeanB = LinePT[OffsetR * 3 + 0];
                    MeanG = LinePT[OffsetR * 3 + 1];
                    MeanR = LinePT[OffsetR * 3 + 2];
                }

                if (Min > LinePSB[OffsetL])
                {
                    Min = LinePSB[OffsetL];
                    MeanB = LinePB[OffsetL * 3 + 0];
                    MeanG = LinePB[OffsetL * 3 + 1];
                    MeanR = LinePB[OffsetL * 3 + 2];
                }

                if (Min > LinePSB[OffsetR])
                {
                    Min = LinePSB[OffsetR];
                    MeanB = LinePB[OffsetR * 3 + 0];
                    MeanG = LinePB[OffsetR * 3 + 1];
                    MeanR = LinePB[OffsetR * 3 + 2];
                }
                LinePD[0] = MeanB;
                LinePD[1] = MeanG;
                LinePD[2] = MeanR;
                LinePD += 3;
            }
        }
    }
FreeMemory:
    if (Average != NULL)        free(Average);
    if (Deviation != NULL)        free(Deviation);
    if (RowOffset != NULL)        free(RowOffset);
    if (ColOffset != NULL)        free(ColOffset);
    if (Blur != NULL)            free(Blur);
    return Status;
}

　　我們跳過內存分配和均方差等的計算過程，我們來看看灰度版本在寬度方向上是如何根據標準函數來獲取值函數的：

            for (int X = 0; X < Width; X++)
            {
                int OffsetL = RowOffset[X];                                            //    求4個點的均方差的最小值
                int OffsetR = RowOffset[X + Radius + Radius];
                int Min, Mean;
                if (LinePST[OffsetL] < LinePST[OffsetR])
                {
                    Min = LinePST[OffsetL];
                    Mean = LinePT[OffsetL];
                }
                else
                {
                    Min = LinePST[OffsetR];
                    Mean = LinePT[OffsetR];
                }
                if (Min > LinePSB[OffsetL])
                {
                    Min = LinePSB[OffsetL];
                    Mean = LinePB[OffsetL];
                }
                if (Min > LinePSB[OffsetR])
                {
                    Min = LinePSB[OffsetR];
                    Mean = LinePB[OffsetR];
                }
                LinePD[X] = Mean;
            }

　　很簡單，就是對四個取樣點，獲取方差的最小值，然後同步獲取對應的均值的值，作為結果值。

本質上講，這個也是個邊緣保留濾波器，雖然其保邊效果不怎麽樣，其實Kuwahara主要是作為一個藝術化的濾鏡存在，在開源的GPUImage裏也有Kuwahara濾波器的實現，參見GPUImageKuwaharaFilter.h文件。

　　我也用別人用的一個美女圖，貼個效果吧。

　　　　　　　　　 　　　原圖　　　　　　　　　　　　　　　　　　　　　　　　　　　 　半徑1　　　　　　　　　　　　　　　　　　　　　　　　 半徑5

　 註意，對於彩色圖像，我們不易將他們分通道來計算，因為分通道時，同一個像素各通道對應的方差最小值不一定在同一位置，此時就會在最終的結果圖像上出現色斑，這不是我們所希望的，在作者的文章中也有類似的說法，如下所示，此時我們可取明度值作為各通道的統一判斷依據。

Obviously the Normal filter can‘t be used for color images by applying the filter to each RGB channel separately and then using the three resulting channels to compose the image. The main problem with this is that the sub regions will have different variancesfor each of the channels. For example, a region with the lowest variance in the red channel might have the highest variance in the green channel. This once again causes abiguity which would result in the color of the central pixel to be determined by several regions, which might also result in blurrier edges.

二、MLV濾波器

　　MLV濾波器也是一個典型的value-and-criterion structure結構濾波器，其參考論文見： A morphology-based filter structure for edge-enhanceing smoothing，和Kuwahara不同的是，MLV濾波器並不只取周邊四領域的均方差作為標準，而是周邊所有包含了中心點店像素的領域的均方差作為判斷依據，取均方差最小那個位置的均值來替代模板中心的值。那這個時候計算的耗時因素就必須納入考慮範圍，因為隨著半徑的增大，這個計算量急劇上升。

　　此時我們想到了前面我的博文SSE圖像算法優化系列七：基於SSE實現的極速的矩形核腐蝕和膨脹（最大值和最小值）算法 中提到了和半徑無關的最值算法的優化，那裏巧妙的運用了一些特性，使得算法在半徑大之後有可能計算時間還有所下降。

　　但是我們註意到這裏的MLV濾波器並不是簡單的最值計算，其還帶了一個尾巴，在計算最值得時候還需要保留另外一個參數同步，這就有點像C的sort函數功能，比如我要對一個結構體數組進行排序，而排序的規則是根據結構體中某一個成員的值來決定的，sort函數可以輕松搞定這個功能，如果我們用普通的C語言實現上面博文的算法，也是可以容易搞定這個的，但是我這裏需要進行SSE處理，那我們先來看看源博文這一塊的核心實現語句。

    __m128i Max1 = _mm_set1_epi8(255), Max2 = _mm_set1_epi8(255);                //    這樣寫能減少一次內存加載
    for (int Y = StartY; Y < EndY; Y++)
    {
        Max1 = _mm_min_epu8(_mm_loadu_si128((__m128i *)(LinePS + 0)), Max1);    //    一條AVX指令
        Max2 = _mm_min_epu8(_mm_loadu_si128((__m128i *)(LinePS + 16)), Max2);
        _mm_store_si128((__m128i *)(LinePD + 0), Max1);
        _mm_store_si128((__m128i *)(LinePD + 16), Max2);
        LinePS += Stride;
        LinePD += 32;
    }

　　很簡單，就是_mm_min_epu8的一個簡單調用，這裏沒有任何的比較函數。

如果我們在計算最值得時候還需要利用最值得特性附帶上另外一個參量，很明顯，如果直接使用上述過程，我們無法做到這一點。必須想辦法帶上比較特征。

　　實際上，min系列的比較函數，我們也可以使用cmp系列的SSE函數來實現，比如要實現兩個32位數的SSE函數的比較，除了使用_mm_min_epi32外，還可以使用下面的代碼：

    __m128i S1 = _mm_set_epi32(10, 20, 30, 50);
    __m128i S2 = _mm_set_epi32(15, 8, 34, 234);

    __m128i Min1 = _mm_min_epi32(S1, S2);

    __m128i Cmp = _mm_cmplt_epi32(S1, S2);
    __m128i Min2 = _mm_blendv_epi8(S2, S1, Cmp);

    for (int Y = 0; Y < 4; Y++)
    {
        printf("%d  \t", Min1.m128i_i32[Y]);
    }
    printf("\n");
    for (int Y = 0; Y < 4; Y++)
    {
        printf("%d  \t", Min2.m128i_i32[Y]);
    }

　　結果如下，完全一樣：

　　註意此時我們在過程中獲得了一些標誌信息，比如上述代碼中得Cmp變量，這個非常有用，有了它，我們就相當於有了一枚旗幟，在隊伍A裏有人舉起了一面旗幟，我們只要在隊伍B裏對應的位置找到對應的人就可以了。這個功能同樣是可以借助於_mm_blendv_epi8這個來實現的。

　　具體到本例，一般方差（不是均方差，此例為未除以取樣總數的值，即(v0-m)*(v0-m)+(v1-m)*(v1-m)+...+(vn-m)*(vn-m)的值，其中m表示v0到vn的平均值，中間用浮點保存），可以用整形表示（在計算最後舍入），而平均值直接使用unsigned char表示，如上IM_KuwaharaFilter函數所示。 此時，按照最值那篇文章的例子，我們可以一次性比較四個SSE變量，這樣產生了4個比較標誌位，而此時我們需要附帶的信息是4*4個字節類型的數據，此時還需將4個比較標誌位的SSE變量合並成一個SSE變量，以方便_mm_blendv_epi8使用，這個時候我們就可以借用packs系列的函數了，具體如下所示:

    __m128i Min0 = _mm_set1_epi32(IM_INT_MAX), Min1 = _mm_set1_epi32(IM_INT_MAX);        //    同時處理16個方差數據
    __m128i Min2 = _mm_set1_epi32(IM_INT_MAX), Min3 = _mm_set1_epi32(IM_INT_MAX);
    __m128i Avg = _mm_setzero_si128(), AvgB = _mm_setzero_si128();
    __m128i AvgG = _mm_setzero_si128(), AvgR = _mm_setzero_si128();
    for (int Y = StartY; Y < EndY; Y++)
    {
        __m128i Src0 = _mm_loadu_si128((__m128i *)(LinePS + 0));
        __m128i Src1 = _mm_loadu_si128((__m128i *)(LinePS + 4));
        __m128i Src2 = _mm_loadu_si128((__m128i *)(LinePS + 8));
        __m128i Src3 = _mm_loadu_si128((__m128i *)(LinePS + 12));
        __m128i Cmp0 = _mm_cmplt_epi32(Src0, Min0);
        __m128i Cmp1 = _mm_cmplt_epi32(Src1, Min1);
        __m128i Cmp2 = _mm_cmplt_epi32(Src2, Min2);
        __m128i Cmp3 = _mm_cmplt_epi32(Src3, Min3);

        Min0 = _mm_blendv_epi8(Min0, Src0, Cmp0);            //    不使用_mm_min_epi32，主要是為了得到標識符
        Min1 = _mm_blendv_epi8(Min1, Src1, Cmp1);
        Min2 = _mm_blendv_epi8(Min2, Src2, Cmp2);
        Min3 = _mm_blendv_epi8(Min3, Src3, Cmp3);
                                                                    
        _mm_storeu_si128((__m128i *)(LinePD + 0), Min0);
        _mm_storeu_si128((__m128i *)(LinePD + 4), Min1);
        _mm_storeu_si128((__m128i *)(LinePD + 8), Min2);
        _mm_storeu_si128((__m128i *)(LinePD + 12), Min3);

        if (Channel == 1)                                //    註意要使用_mm_packs_epi16而不是_mm_packus_epi16
        {
            Avg = _mm_blendv_epi8(Avg, _mm_loadu_si128((__m128i *)LinePA), _mm_packs_epi16(_mm_packs_epi32(Cmp0, Cmp1), _mm_packs_epi32(Cmp2, Cmp3)));
            _mm_storeu_si128((__m128i *)LinePM, Avg);            //    最小方差對應的均值
        }
        else
        {
            __m128i Mask0 = _mm_or_si128(_mm_shuffle_epi8(Cmp0, _mm_setr_epi8(0, 0, 0, 4, 4, 4, 8, 8, 8, 12, 12, 12, -1, -1, -1, -1)), _mm_shuffle_epi8(Cmp1, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 4)));
            __m128i Mask1 = _mm_or_si128(_mm_shuffle_epi8(Cmp1, _mm_setr_epi8(4, 4, 8, 8, 8, 12, 12, 12, -1, -1, -1, -1, -1, -1, -1, -1)), _mm_shuffle_epi8(Cmp2, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 4, 4, 4, 8, 8)));
            __m128i Mask2 = _mm_or_si128(_mm_shuffle_epi8(Cmp2, _mm_setr_epi8(8, 12, 12, 12, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1)), _mm_shuffle_epi8(Cmp3, _mm_setr_epi8(-1, -1, -1, -1, 0, 0, 0, 4, 4, 4, 8, 8, 8, 12, 12, 12)));
            AvgB = _mm_blendv_epi8(AvgB, _mm_loadu_si128((__m128i *)LinePA), Mask0);
            AvgG = _mm_blendv_epi8(AvgG, _mm_loadu_si128((__m128i *)(LinePA + 16)), Mask1);
            AvgR = _mm_blendv_epi8(AvgR, _mm_loadu_si128((__m128i *)(LinePA + 32)), Mask2);
            _mm_storeu_si128((__m128i *)LinePM, AvgB);
            _mm_storeu_si128((__m128i *)(LinePM + 16), AvgG);
            _mm_storeu_si128((__m128i *)(LinePM + 32), AvgR);
        }
        LinePS += Width;
        LinePA += Stride;
        LinePD += 16;
        LinePM += 16 * Channel;
    }

　　註意在上述代碼中我們使用packs而非packus，請讀者自行理解這是為什麽。

　　對於彩色圖像，根據上面的對Kuwahara濾波器的描述，也是需要計算明度的，此時，一個標誌位對應的就是就是3個分量，此時就不能使用packs之類的函數，而可以直接使用上述的_mm_shuffle_epi8類的函數可以輕松搞定。

　　當然具體的過程還涉及到很多其他東西，比如int類型矩陣的轉置，內存的有效利用等等，這裏不予以多述。

　　三、MCV濾波器

　　同MLV濾波器稍微有點不同的時，MCV濾波器的評判標準有稍作改動，其參考文件見：An edge-enhancing nonlinear filter for reducing multiplicative noise，但是核心的實現沒有太大的區別。

　　按照論文的說法，MCV濾波器可以有效地去除乘性噪音，而MLV可以去除加性噪音，如下所示：

The MCV filter therefore uses the sample mean for the value function, the coefficient of variation as the criterion function, and the minimum as the selection function. It is a value-and-criterion filter specifically designed to reduce multiplicative noise.

　　以及

The MCV filter is closely related to a value-and-criterion filter designed to remove additive noise known as the Mean of Least Variance (MLV) filter [7–9]. The MLV filter uses the sample variance as its criterion function rather than the coefficient of variation. In images corrupted by additive noise, the variance is theoretically minimum in structuring elements where the signal is constant.

　　MLV和MCV能產生比Kuwahara更為平滑的圖像，在SAR圖像以及醫療圖像上應該有一定用武之地。而且個人感覺在小半徑時其產生的結果對圖像的邊界分割、邊緣查找等也有一定的幫助。

　　　　　　　　　　　　　　　　 　　原圖　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 　Kuwahara

　　　　　　　　　　　　　　　　　　 　　MLV　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　MCV

　　上述操作半徑都為5。

　　處理後，似乎邊緣的分界線更為明顯。

　　當然我貼的圖形其實都不是論文本身強調的應用場景，只是做個例證而已。

　　關於這種類似的濾波器，還可以使用圓形半徑，或者是橢圓半徑，也有一些人對這方面的需求做了研究，比如https://blog.csdn.net/qq_16013649/article/details/78744910。

　　針對1024*768的彩色圖像，此類算法優化後大約耗時20ms，灰度圖約10ms，速度也還算比較快了，測試見附件的SSE_Optimization_Demo的stylize菜單。

Demo下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

SSE圖像算法優化系列二十三: 基於value-and-criterion structure 系列濾波器（如Kuwahara，MLV，MCV濾波器）的優化。