當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

加权中值滤波

發布時間：2025/3/17 编程问答 14 豆豆

生活随笔收集整理的這篇文章主要介紹了加权中值滤波小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

應作者要求，轉載自：https://www.cnblogs.com/Imageshop/p/9934670.html

《100+ Times FasterWeighted Median Filter (WMF)》

　　這篇文章的官網地址是：http://www.cse.cuhk.edu.hk/~leojia/projects/fastwmedian/，其中主要作者Jiaya Jia教授的官網地址是：http://jiaya.me/，根據Jiaya Jia的說法，這個算法很快將被OpenCv所收錄，到時候OpenCv的大神應該對他還有所改進吧。

　　在百度上搜索加權中值模糊，似乎只有一篇博客對這個文章進行了簡單的描述，詳見：https://blog.csdn.net/streamchuanxi/article/details/79573302?utm_source=blogxgwz9。

　　由于作者只給出了最后的優化實現代碼，而論文中還提出了各種中間過程的時間，因此本文以實現和驗證論文中有關說法為主，涉及到的理論知識比較膚淺，一般是一筆而過。

　　根據論文中得說法，所謂的加權中值濾波，也是一種非線性的圖像平滑技術，他取一個局部窗口內所有像素的加權中值來代替局部窗口的中心點的值。用較為數學的方法表示如下：

　　在圖像I中的像素p，我們考慮以p為中心，半徑為R的局部窗口，不同于普通的中值模糊，對于屬于內每一個像素q，都有一個基于對應的特征圖像的相似度的權重系數wpq，如下式所示：

? ? ? ? ? ? ? ? ? ? ? ? ? ?

　　f(p)和f(q)是像素p和q在對應的特征圖中得特征值。g是一個權重函數，最常用的即為高斯函數，反應了像素p和q的相似程度。

　　我們用I(q)表示像素點q的像素值，在窗口內的像素總數量用n表示，則n=(2r+1)*(2r+1)，那么窗口內像素值和權重值構成一個對序列，即，對這個序列按照I(q)的值進行排序。排序后，我們依次累加權重值，直到累加的權重大于等于所有權重值的一半時停止，此時對應的I(q)即作為本局部窗口中心點的新的像素值。

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

　　很明顯，上面的過程要比標準的中值模糊復雜一些，在處理時多了特征圖和權重函數項，而標準的中值模糊我們可以認為是加權中值模糊的特例，即所有局部窗口的權重都為1或者說相等。

　　在這里，特征圖可以直接是源圖像，也可以是其他的一些特征，比如原圖像的邊緣檢測結果、局部均方差、局部熵或者其他的更為高級的特征。

　　按照這個定義，我們給出一段針對灰度數據的Brute-force處理代碼：

int __cdecl ComparisonFunction(const void *X, const void *Y) // 一定要用__cdecl這個標識符 {Value_Weight VWX = *(Value_Weight *)X;Value_Weight VWY = *(Value_Weight *)Y;if (VWX.Value < VWY.Value)return -1;else if (VWX.Value > VWY.Value)return +1;elsereturn 0; }// 加權中值模糊，直接按照算法的定義實現。 // Input - 輸入圖像，灰度圖，LevelV = 256級 // FeatureMap - 特征圖像，灰度圖，LevelF = 256級 // Weight - 特征的權重矩陣，大小是LevelF * LevelF // Output - 輸出圖像，不能和Input為同一個數據。int IM_WeightedMedianBlur_00(unsigned char *Input, unsigned char *FeatureMap, float *Weight, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((FeatureMap == NULL) || (Weight == NULL)) return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1)) 　　　　　　　　　　return IM_STATUS_NOTSUPPORTED;const int LevelV = 256; // Value 可能出現的不同數量const int LevelF = 256; // Feature 可能出現的不同數量Value_Weight *VW = (Value_Weight *)malloc((2 * Radius + 1) * (2 * Radius + 1) * sizeof(Value_Weight)); // 值和特征序列對內存if (VW == NULL) return IM_STATUS_OK;for (int Y = 0; Y < Height; Y++){unsigned char *LinePF = FeatureMap + Y * Stride;unsigned char *LinePD = Output + Y * Stride;for (int X = 0; X < Width; X++){int CF_Index = LinePF[X] * LevelF;int PixelAmount = 0;float SumW = 0;for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Index = J * Stride;for (int I = IM_Max(X - Radius, 0); I <= IM_Min(X + Radius, Width - 1); I++) // 注意越界{int Value = Input[Index + I]; // 值int Feature = FeatureMap[Index + I]; // 特征float CurWeight = Weight[CF_Index + Feature]; // 對應的權重VW[PixelAmount].Value = Value;VW[PixelAmount].Weight = CurWeight; // 保存數據SumW += CurWeight; // 計算累加數據PixelAmount++; // 有效的數據量 }}float HalfSumW = SumW * 0.5f; // 一半的權重SumW = 0;qsort(VW, PixelAmount, sizeof VW[0], &ComparisonFunction); // 調用系統的qsort按照Value的值從小到大排序，注意qsort的結果仍然保存在第一個參數中for (int I = 0; I < PixelAmount; I++) // 計算中值{SumW += VW[I].Weight;if (SumW >= HalfSumW){LinePD[X] = VW[I].Value;break;}}}}free(VW);return IM_STATUS_OK; }

　　很明顯，這個函數的時間復雜度是o(radius * radius)，空間復雜度到時很小。

　　我們在一臺 I5，3.3GHZ的機器上進行了測試，上述代碼處理一副1000*1000像素的灰度圖，半徑為10（窗口大小21*21）時，處理時間約為27s，論文里給的Cpu和我的差不多，給出的處理one - metalpixel的RGB圖用時90.7s，考慮到RGB的通道的數據量以及一些其他的處理，應該說論文如實匯報了測試數據。

　　那么從代碼優化上面講，上面代碼雖然還有優化的地方，但是都是小打小鬧了。使用VS的性能分析器，可以大概獲得如下的結果：

? ? ? ?

　　可見核心代碼基本都用于排序了，使用更快的排序有助于進一步提高速度。

　　針對這個情況，論文的作者從多方面提出了改進措施，主要有三個方面，我們簡單的重復下。

　　一、聯合直方圖(Joint Histgram）

　　直方圖優化在很多算法中都有應用，比如標準的中值濾波，現在看到的最快的實現方式還是基于直方圖的，詳見：任意半徑中值濾波（擴展至百分比濾波器）O(1)時間復雜度算法的原理、實現及效果，但是在加權中值濾波中，傳統的一維直方圖已經無法應用，因為這個算法不僅涉及到原圖的像素值，還和另外一幅特征圖有關，因此，文中提出了聯合直方圖，也是一種二維直方圖。

　　如果圖像中的像素最多有LevelV個不同值，其對應的特征最多有LevelF個不同的值，那么我們定義一個寬和高分別為LevelV * LevelF大小的直方圖。對于某一個窗口，統計其內部的(2r+1)*(2r+1)個像素和特征對的直方圖數據，即如果某個點的像素值為V，對應的特征值為F，則相應位置的直方圖數據加1。

　　如果我們統計出這個二維的直方圖數據后，由于中心點的特征值是固定的，因此，對于直方圖的每一個LevelF值，權重是一定的了，我們只需計算出直方圖內每一個Value值所對應所有的Feature的權重后，就可方便的統計出中值所在的位置了。

　　那么如果每個像素點都進行領域直方圖的計算，這個的工作量也是蠻大的，同一維直方圖的優化思路一樣，在進行逐像素行處理的時候，對直方圖數據可以進行逐步的更新，去除掉移走的那一列的直方圖信息，在加入即將進入那一列數據，而中間重疊部分則不需要調整。

　　按照論文中的Joint Histgram的布局，即行方向大小為LevelV，列方向大小為LevelF，編制Joint Histgram實現的加權中值算法代碼如下所示：

// 加權中值模糊，基于論文中圖示的內存布局設置的Joint Histgram。 int IM_WeightedMedianBlur_01(unsigned char *Input, unsigned char *FeatureMap, float *Weight, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((FeatureMap == NULL) || (Weight == NULL)) return IM_STATUS_NULLREFRENCE; if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1) && (Channel != 3)) return IM_STATUS_NOTSUPPORTED;int Status = IM_STATUS_OK;const int LevelV = 256; // Value 可能出現的不同數量const int LevelF = 256; // Feature 可能出現的不同數量int *Histgram = (int *)malloc(LevelF * LevelV * sizeof(int));float *Sum = (float *)malloc(LevelV * sizeof(float));if ((Histgram == NULL) || (Sum == NULL)){Status = IM_STATUS_OUTOFMEMORY;goto FreeMemory;}for (int Y = 0; Y < Height; Y++){unsigned char *LinePF = FeatureMap + Y * Stride;unsigned char *LinePD = Output + Y * Stride;memset(Histgram, 0, LevelF * LevelV * sizeof(int));for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){for (int I = IM_Max(0 - Radius, 0); I <= IM_Min(0 + Radius, Width - 1); I++){int Value = Input[J * Stride + I];int Feature = FeatureMap[J * Stride + I]; // 統計二維直方圖Histgram[Feature * LevelV + Value]++;}}for (int X = 0; X < Width; X++){int Feature = LinePF[X];float SumW = 0, HalfSumW = 0;;for (int I = 0; I < LevelV; I++){float Cum = 0;for (int J = 0; J < LevelF; J++) // 計算每個Value列針對的不同的Feature的權重的累計值{Cum += Histgram[J * LevelV + I] * Weight[J * LevelF + Feature];}Sum[I] = Cum;SumW += Cum;}HalfSumW = SumW / 2;SumW = 0;for (int I = 0; I < LevelV; I++){SumW += Sum[I];if (SumW >= HalfSumW) // 計算中值{LinePD[X] = I;break;}}if ((X - Radius) >= 0) // 移出的那一列的直方圖{for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X - Radius];int Feature = FeatureMap[J * Stride + X - Radius];Histgram[Feature * LevelV + Value]--;}}if ((X + Radius + 1) <= Width - 1) // 移入的那一列的直方圖{for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X + Radius + 1];int Feature = FeatureMap[J * Stride + X + Radius + 1];Histgram[Feature * LevelV + Value]++;}}}} FreeMemory:if (Histgram != NULL) free(Histgram);if (Sum != NULL) free(Sum);return Status; }

　　編譯后測試，同樣是21*21的窗口，one - metalpixel的灰度圖像計算用時多達108s，比直接實現慢很多了。

　　分析原因，核心就是在中值的查找上，由于我們采用的內存布局方式，導致計算每個Value對應的權重累加存在的大量的Cache miss現象，即下面這條語句：

for (int J = 0; J < LevelF; J++) // 計算每個Value列針對的不同的Feature的權重的累計值 {Cum += Histgram[J * LevelV + I] * Weight[J * LevelF + Feature]; }

　　我們換種Joint Histgram的布局，即行方向大小為LevelF，列方向大小為LevelV，此時的代碼如下：

// 加權中值模糊，修改內存布局設置的Joint Histgram。 int IM_WeightedMedianBlur_02(unsigned char *Input, unsigned char *FeatureMap, float *Weight, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((FeatureMap == NULL) || (Weight == NULL)) return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1) && (Channel != 3)) return IM_STATUS_NOTSUPPORTED;int Status = IM_STATUS_OK;const int LevelV = 256; // Value 可能出現的不同數量const int LevelF = 256; // Feature 可能出現的不同數量int *Histgram = (int *)malloc(LevelF * LevelV * sizeof(int));float *Sum = (float *)malloc(LevelV * sizeof(float));if ((Histgram == NULL) || (Sum == NULL)){Status = IM_STATUS_OUTOFMEMORY;goto FreeMemory;}for (int Y = 0; Y < Height; Y++){unsigned char *LinePF = FeatureMap + Y * Stride;unsigned char *LinePD = Output + Y * Stride;memset(Histgram, 0, LevelF * LevelV * sizeof(int));for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Index = J * Stride;for (int I = IM_Max(0 - Radius, 0); I <= IM_Min(0 + Radius, Width - 1); I++){int Value = Input[J * Stride + I];int Feature = FeatureMap[J * Stride + I];Histgram[Value * LevelF + Feature]++; // 注意索引的方式的不同}}for (int X = 0; X < Width; X++){int IndexF = LinePF[X] * LevelF;float SumW = 0, HalfSumW = 0;;for (int I = 0; I < LevelV; I++){float Cum = 0;int Index = I * LevelF;for (int J = 0; J < LevelF; J++) // 核心就這里不同{Cum += Histgram[Index + J] * Weight[IndexF + J];}Sum[I] = Cum;SumW += Cum;}HalfSumW = SumW / 2;SumW = 0;for (int I = 0; I < LevelV; I++){SumW += Sum[I];if (SumW >= HalfSumW){LinePD[X] = I;break;}}if ((X - Radius) >= 0){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X - Radius];int Feature = FeatureMap[J * Stride + X - Radius];Histgram[Value * LevelF + Feature]--;}}if ((X + Radius + 1) <= Width - 1){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X + Radius + 1];int Feature = FeatureMap[J * Stride + X + Radius + 1];Histgram[Value * LevelF + Feature]++;}}}} FreeMemory:if (Histgram != NULL) free(Histgram);if (Sum != NULL) free(Sum);return Status; }

　　修改后，同樣的測試條件和圖片，速度提升到了17s，僅僅是更改了一個內存布局而已，原論文的圖沒有采用這種布局方式，也許只是為了表達算法清晰而已。

　　和原論文比較，原論文的joint histgram時間要比直接實現慢（156.9s vs 90.7s)，而我這里的一個版本比brute force的快，一個比brute force的慢，因此，不清楚作者在比較時采用了何種編碼方式，但是這都不重要，因為他們的區別都還在一個數量級上。

? ? ? ?由于直方圖大小是固定的，因此，前面的中值查找的時間復雜度是固定的，而后續的直方圖更新則是o(r)的，但是注意到由于LevelV和 LevelF通常都是比較大的常數（一般為256），因此實際上，中值查找這一塊的耗時占了絕對的比例。

　二、快速中值追蹤

　　尋找中值的過程實際上可以看成一個追求平衡的過程，假定當前搜索到的位置是V，位于V左側所有相關值的和是Wl，位于V右側所有相關值得和是Wr，則中值的尋找可以認為是下式：

? ? ? ? ? ? ? ? ? ? ? ? ??

　　后面的約束條件可以理解為第一次出現Wl大于Wr前。

? ? ? ?如果我們之前已經尋找到了像素P處的中值，那么由于像素的連續性，像素P+1處的中值一般不會和P處的中值差異太大，大量的統計數據表明他們的差異基本在8個像素值之類（256色階圖），那么這個思想其實和任意半徑中值濾波（擴展至百分比濾波器）O(1)時間復雜度算法的原理、實現及效果中講到的是一致的。這種特性，我們也可以將他運用到加權中值濾波中。

　　考慮到加權中值濾波中聯合直方圖的特殊性，我們需要維護一張平衡表，論文中叫做Balance Counting Box(BCB)，這一塊的解釋比較拗口也比較晦澀，大家需要仔細的看論文和我下面提供的JointHist+MedianTracking代碼。

// 加權中值模糊, Joint + MT int IM_WeightedMedianBlur_03(unsigned char *Input, unsigned char *FeatureMap, float *Weight, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((FeatureMap == NULL) || (Weight == NULL)) return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1) && (Channel != 3)) return IM_STATUS_NOTSUPPORTED;int Status = IM_STATUS_OK;const int LevelV = 256; // Value 可能出現的不同數量const int LevelF = 256; // Feature 可能出現的不同數量int *Histgram = (int *)malloc(LevelF * LevelV * sizeof(int));int *BCB = (int *)malloc(LevelF * sizeof(int));if ((Histgram == NULL) || (BCB == NULL)){Status = IM_STATUS_OK;return IM_STATUS_OUTOFMEMORY;}for (int Y = 0; Y < Height; Y++){unsigned char *LinePF = FeatureMap + Y * Stride;unsigned char *LinePD = Output + Y * Stride;memset(Histgram, 0, LevelF * LevelV * sizeof(int)); // 全部賦值為0memset(BCB, 0, LevelF * sizeof(int));int CutPoint = -1;for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Index = J * Stride;for (int I = IM_Max(0 - Radius, 0); I <= IM_Min(0 + Radius, Width - 1); I++){int Value = Input[J * Stride + I];int Feature = FeatureMap[J * Stride + I];Histgram[Value * LevelF + Feature]++; // 計算每行第一個點的二維直方圖，直方圖的水平方向為Feature坐標，垂直方向為Value坐標 BCB[Feature]--; // 此時的CutPoint初始化為-1，所以+方向的數據為0，所有的都在-方向 }}for (int X = 0; X < Width; X++){float BalanceWeight = 0;int IndexF = LinePF[X] * LevelF; // 中心點P的Value所對應的那一行Feature權重起始索引for (int I = 0; I < LevelF; I++) // BCB[I]中保存的是以CutPoint為分界線，Feature為I時，分界線左側的所有Value[0-CutPoint]值的數量和分界線右側所有的Value（CutPoint， LevelV - 1]值數量的差異{BalanceWeight += BCB[I] * Weight[IndexF + I]; // 因為Feature為固定值時，如果中心點固定，那么不管與Feature對應的Value值時多少，Weight就是定值了。}if (BalanceWeight < 0) // 第一個點的BalanceWeight必然小于0{for (; BalanceWeight < 0 && CutPoint != LevelV - 1; CutPoint++){int IndexH = (CutPoint + 1) * LevelF; // 新的直方圖的位置float CurWeight = 0;for (int I = 0; I < LevelF; I++){CurWeight += 2 * Histgram[IndexH + I] * Weight[IndexF + I]; // 左側加右側同時減，所以是2倍BCB[I] += Histgram[IndexH + I] * 2; // 數量是同樣的道理}BalanceWeight += CurWeight;}}else if (BalanceWeight > 0) // 如果平衡值大于0，則向左移動中間值{for (; BalanceWeight > 0 && CutPoint != 0; CutPoint--){int IndexH = CutPoint * LevelF;float CurWeight = 0;for (int I = 0; I < LevelF; I++){CurWeight += 2 * Histgram[IndexH + I] * Weight[IndexF + I];BCB[I] -= Histgram[IndexH + I] * 2;}BalanceWeight -= CurWeight;}}LinePD[X] = CutPoint;if ((X - Radius) >= 0){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++) // 即將移出的那一列數據{int Value = Input[J * Stride + X - Radius];int Feature = FeatureMap[J * Stride + X - Radius];Histgram[Value * LevelF + Feature]--;if (Value <= CutPoint) // 如果移出的那個值小于當前的中值BCB[Feature]--;elseBCB[Feature]++;}}if ((X + Radius + 1) <= Width - 1){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X + Radius + 1];int Feature = FeatureMap[J * Stride + X + Radius + 1];Histgram[Value * LevelF + Feature]++;if (Value <= CutPoint) // 如果移出的那個值小于當前的中值BCB[Feature]++;elseBCB[Feature]--;}}}}free(Histgram);free(BCB); }

　　代碼也很簡潔，主要是增加了一個BCB列表的維護，編譯后測試，同樣是21*21的窗口，one - metalpixel的灰度圖像計算用420ms左右，比Brute-force版本的27s大約快了64倍，這個和論文的時間比例基本差不多(（156.9+0.4）/(2.2+0.5)=58)。提速也是相當的可觀，而且算法速度和半徑不是特別敏感，畢竟更新直方圖的計算量在這里占的比例其實已經不多了。

　　三、Necklace Table

? ? 那么論文最后還提出了另外的進一步加速的方案，這是基于以下觀察到的事實，即在直方圖的數據中，存在大量的0值，這些值的計算其實對算法本身是沒有任何作用的，但是占用了大量的計算時間。

? ??

　　比如上圖是某個圖像局部窗口的聯合直方圖和BCB值，在聯合直方圖中大部分區域都是0值對應的黑色，在BCB中大部分情況也是0值。

? ? ? ?因此，作者構建了一個叫做Necklace Table的數據結構，這個數據結構可以方便快捷的記錄下一個和上一個非0元素的位置，從而能有效的訪問到那些真正有計算價值的部位，以及簡單的刪除和增加節點的功能，具體的實現細節詳見論文或下面的JointHistgram + Necklace Table代碼。

int IM_WeightedMedianBlur_04(unsigned char *Input, unsigned char *FeatureMap, float *Weight, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((FeatureMap == NULL) || (Weight == NULL)) return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1) && (Channel != 3)) return IM_STATUS_NOTSUPPORTED;int Status = IM_STATUS_OK;const int LevelV = 256; // Value 可能出現的不同數量const int LevelF = 256; // Feature 可能出現的不同數量 const int LevelV = 256;int *Histgram = (int *)malloc(LevelF * LevelV * sizeof(int));int *ForwardH = (int *)malloc(LevelF * LevelV * sizeof(int)); // forward link for necklace tableint *BackWordH = (int *)malloc(LevelF * LevelV * sizeof(int)); // forward link for necklace tablefloat *Sum = (float *)malloc(LevelV * sizeof(float));if ((Histgram == NULL) || (ForwardH == NULL) || (BackWordH == NULL) || (Sum == NULL)){Status = IM_STATUS_OK;goto FreeMemory;}memset(ForwardH, 0, LevelF * LevelV * sizeof(int));memset(BackWordH, 0, LevelF * LevelV * sizeof(int));for (int Y = 0; Y < Height; Y++){unsigned char *LinePF = FeatureMap + Y * Stride;unsigned char *LinePD = Output + Y * Stride;memset(Histgram, 0, LevelF * LevelV * sizeof(int));for (int X = 0; X < LevelV; X++){ForwardH[X * LevelF] = 0; // 其實每一個Feature對應一個完整的Necklace Table，需要把第一個元素置為0BackWordH[X * LevelF] = 0;}for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++) // 第一個元素{int Index = J * Stride;for (int I = IM_Max(0 - Radius, 0); I <= IM_Min(0 + Radius, Width - 1); I++){int Value = Input[Index + I];int Feature = FeatureMap[Index + I];int Index = Value * LevelF;if (Histgram[Index + Feature] == 0 && Feature != 0) // 直方圖數據如果還是0并且FMap值不為0{int T = ForwardH[Index];ForwardH[Index] = Feature;ForwardH[Index + Feature] = T;BackWordH[Index + T] = Feature;BackWordH[Index + Feature] = 0;}Histgram[Index + Feature]++;}}for (int X = 0; X < Width; X++){int IndexF = LinePF[X] * LevelF;float SumW = 0, HalfSumW = 0;;for (int I = 0; I < LevelV; I++){float Cum = 0;int Index = I * LevelF;int J = 0;do{Cum += Histgram[Index + J] * Weight[IndexF + J]; // 跳過那些非0的元素J = ForwardH[Index + J];} while (J != 0);Sum[I] = Cum; // 計算每一個Value對應的所有Featrue的權重累計和SumW += Cum;}HalfSumW = SumW / 2;SumW = 0;for (int I = 0; I < LevelV; I++){SumW += Sum[I];if (SumW >= HalfSumW){LinePD[X] = I;break;}}if ((X - Radius) >= 0){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X - Radius];int Feature = FeatureMap[J * Stride + X - Radius];int Index = Value * LevelF;Histgram[Index + Feature]--;if (Histgram[Index + Feature] == 0 && Feature != 0){int T1 = BackWordH[Index + Feature];int T2 = ForwardH[Index + Feature];ForwardH[Index + T1] = T2;BackWordH[Index + T2] = T1;}}}if ((X + Radius + 1) <= Width - 1){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X + Radius + 1];int Feature = FeatureMap[J * Stride + X + Radius + 1];int Index = Value * LevelF;if (Histgram[Index + Feature] == 0 && Feature != 0) // 直方圖數據如果還是0并且FMap值不為0{int T = ForwardH[Index];ForwardH[Index] = Feature;ForwardH[Index + Feature] = T;BackWordH[Index + T] = Feature;BackWordH[Index + Feature] = 0;}Histgram[Index + Feature]++;}}}} FreeMemory:if (Histgram != NULL) free(Histgram);if (ForwardH != NULL) free(ForwardH);if (BackWordH != NULL) free(BackWordH);if (Sum != NULL) free(Sum);return Status; }

? ? ? 代碼量不大，編譯后測試，同樣是21*21的窗口，one - metalpixel的灰度圖像計算用1200ms左右，比Brute-force版本的27s大約快了22倍，由于這個算法和圖像內容是由一定關系的，因此，和論文提供的數據直接比較的意義不大。

? ? ?四、最終的結合體

　很自然的，我們想到要把Median Tracking 和 Necklace Table聯合在一起，來進一步的提高速度，這個時候可以對Joint Histgram即BCB都使用?Necklace Table來記錄非零元素，于是產生了以下的結合代碼：

int IM_WeightedMedianBlur_05(unsigned char *Input, unsigned char *FeatureMap, float *Weight, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((FeatureMap == NULL) || (Weight == NULL)) return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1) && (Channel != 3) && (Channel != 4)) return IM_STATUS_NOTSUPPORTED;int Status = IM_STATUS_OK;const int LevelV = 256;const int LevelF = 256;int *Histgram = (int *)malloc(LevelF * LevelV * sizeof(int));int *BCB = (int *)malloc(LevelF * sizeof(int));int *ForwardH = (int *)malloc(LevelF * LevelV * sizeof(int)); // forward link for necklace tableint *BackWordH = (int *)malloc(LevelF * LevelV * sizeof(int)); // forward link for necklace tableint *ForwardBCB = (int *)malloc(LevelF * sizeof(int)); // forward link for necklace tableint *BackWordBCB = (int *)malloc(LevelF * sizeof(int)); // forward link for necklace tableif ((Histgram == NULL) || (BCB == NULL) || (ForwardH == NULL) || (BackWordH == NULL) || (ForwardBCB == NULL) || (BackWordBCB == NULL)){Status = IM_STATUS_OK;goto FreeMemory;}memset(ForwardH, 0, LevelF * LevelV * sizeof(int));memset(BackWordH, 0, LevelF * LevelV * sizeof(int));memset(ForwardBCB, 0, LevelF * sizeof(int));memset(BackWordBCB, 0, LevelF * sizeof(int));for (int Y = 0; Y < Height; Y++){unsigned char *LinePF = FeatureMap + Y * Stride;unsigned char *LinePD = Output + Y * Stride;memset(Histgram, 0, LevelF * LevelV * sizeof(int)); // 全部賦值為0memset(BCB, 0, LevelF * sizeof(int));for (int X = 0; X < LevelV; X++){ForwardH[X * LevelF] = 0;BackWordH[X * LevelF] = 0;}ForwardBCB[0] = 0;BackWordBCB[0] = 0;int CutPoint = -1;for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Index = J * Stride;for (int I = IM_Max(0 - Radius, 0); I <= IM_Min(0 + Radius, Width - 1); I++){int Value = Input[Index + I];int Feature = FeatureMap[Index + I];int Index = Value * LevelF;if (Histgram[Index + Feature] == 0 && Feature != 0) // 直方圖數據如果還是0并且FMap值不為0{int T = ForwardH[Index];ForwardH[Index] = Feature;ForwardH[Index + Feature] = T;BackWordH[Index + T] = Feature;BackWordH[Index + Feature] = 0;}Histgram[Index + Feature]++; // 計算每行第一個點的二維直方圖，直方圖的水平方向為Feature坐標，垂直方向為Value坐標 UpdateBCB(BCB[Feature], ForwardBCB, BackWordBCB, Feature, -1); // 此時的CutPoint初始化為-1，所以+方向的數據為0，所有的都在-方向 }}for (int X = 0; X < Width; X++){float BalanceWeight = 0;int IndexF = LinePF[X] * LevelF; // 中心點P的Value所對應的那一行Feature權重起始索引int I = 0;do{BalanceWeight += BCB[I] * Weight[IndexF + I]; // 按照當前BCB數據計算平衡值，BCB記錄了相同的FMap值時按照之前的中間值左右兩側像素個數的差異值I = ForwardBCB[I];} while (I != 0);if (BalanceWeight < 0) // 第一個點的BalanceWeight必然小于0{for (; BalanceWeight < 0 && CutPoint != LevelV - 1; CutPoint++){int IndexH = (CutPoint + 1) * LevelF; // 新的直方圖的位置float CurWeight = 0;int I = 0;do{CurWeight += 2 * Histgram[IndexH + I] * Weight[IndexF + I]; // 左側加右側同時減，所以是2倍UpdateBCB(BCB[I], ForwardBCB, BackWordBCB, I, Histgram[IndexH + I] << 1);I = ForwardH[IndexH + I];} while (I != 0);BalanceWeight += CurWeight;}}else if (BalanceWeight > 0) // 如果平衡值大于0，則向左移動中間值{for (; BalanceWeight > 0 && CutPoint != 0; CutPoint--){int IndexH = CutPoint * LevelF;float CurWeight = 0;int I = 0;do{CurWeight += 2 * Histgram[IndexH + I] * Weight[IndexF + I]; // 左側加右側同時減，所以是2倍UpdateBCB(BCB[I], ForwardBCB, BackWordBCB, I, -(Histgram[IndexH + I] << 1));I = ForwardH[IndexH + I];} while (I != 0);BalanceWeight -= CurWeight;}}LinePD[X] = CutPoint;if ((X - Radius) >= 0){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++) // 即將移出的那一列數據{int Value = Input[J * Stride + X - Radius];int Feature = FeatureMap[J * Stride + X - Radius];int Index = Value * LevelF;Histgram[Index + Feature]--;if (Histgram[Index + Feature] == 0 && Feature != 0){int T1 = BackWordH[Index + Feature];int T2 = ForwardH[Index + Feature];ForwardH[Index + T1] = T2;BackWordH[Index + T2] = T1;}UpdateBCB(BCB[Feature], ForwardBCB, BackWordBCB, Feature, -((Value <= CutPoint) << 1) + 1);}}if ((X + Radius + 1) <= Width - 1){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X + Radius + 1];int Feature = FeatureMap[J * Stride + X + Radius + 1];int Index = Value * LevelF;if (Histgram[Index + Feature] == 0 && Feature != 0) // 直方圖數據如果還是0并且FMap值不為0{int T = ForwardH[Index];ForwardH[Index] = Feature;ForwardH[Index + Feature] = T;BackWordH[Index + T] = Feature;BackWordH[Index + Feature] = 0;}UpdateBCB(BCB[Feature], ForwardBCB, BackWordBCB, Feature, ((Value <= CutPoint) << 1) - 1);Histgram[Index + Feature]++;}}}} FreeMemory:if (Histgram != NULL) free(Histgram);if (BCB != NULL) free(BCB);if (ForwardH != NULL) free(ForwardH);if (BackWordH != NULL) free(BackWordH);if (ForwardBCB != NULL) free(ForwardBCB);if (BackWordBCB != NULL) free(BackWordBCB);return Status; }

　　我們滿懷期待的編譯和執行他，結果出來了，同樣是21*21的窗口，one - metalpixel的灰度圖像計算用430ms左右，和Joint + MT的速度差不多，但是論文里給出的數據是Joint + MT + NT要比Joint + MT快3倍左右。這是怎么回事呢。

　　我們仔細檢查論文里，在Implementation Notes節里有這樣的語句：

? ? ? ? ? ? ??Only a single thread is used without involving any SIMD?instructions. Our system is implemented using C++.?

　　第一，他也是用的C++和我一樣，第二，他是單線程，也和我一樣，第三，沒有使用任何SIMD指令，似乎我也沒有使用啊，都一樣，為什么結果比對不一致，難道是大神他們作弊，鑒于他們的成就，我立即撤回我這逆天的想法，一定是其他地方有問題。我們試著反編譯看看。

　　我們定位到Joint + MT的算法的下面一句代碼看看：

　　反編譯結果為：

　　赤裸裸的SIMD指令啊。

　　為什么呢，只是因為VS的編譯器即使在默認情況下的設置中，也會根據當前編譯系統的情況，進行一定的向量化優化，加上現在的PC基本沒有哪一個不能使用SIMD指令的。如下圖所示，為C++默認編譯選項：

? ? ? ? ? ?

　　在啟用增強指令集選項里默認是未設置，但是未設置并不代表不使用，正如上述所言，測試編譯器會根據系統狀況優化編譯。因此，雖然表面上代碼沒有使用SIMD指令，但是實際卻使用了。

　　為了公平起見，我們禁用系統的SIMD優化，此時，可以在增強指令集的選項里選擇“無增強指令/arch:IA32".

? ? ? ? ?

　編譯后，對上述同樣一段代碼進行反編譯，可以看到如下匯編碼：

for (int I = 0; I < LevelF; I++) // BCB[I]中保存的是以CutPoint為分界線，Feature為I時，分界線左側的所有Value[0-CutPoint]值的數量和分界線右側所有的Value（CutPoint， LevelV - 1]值數量的差異{BalanceWeight += BCB[I] * Weight[IndexF + I]; // 因為Feature為固定值時，如果中心點固定，那么不管與Feature對應的Value值時多少，Weight就是定值了。 0F8F1AF5 fild dword ptr [ecx-4] 0F8F1AF8 fmul dword ptr [eax+4] 0F8F1AFB fild dword ptr [ecx-8] 0F8F1AFE fmul dword ptr [eax] 0F8F1B00 faddp st(2),st 0F8F1B02 faddp st(1),st 0F8F1B04 fild dword ptr [ecx] 0F8F1B06 fmul dword ptr [eax+8] 0F8F1B09 faddp st(1),st 0F8F1B0B fild dword ptr [ecx+4] 0F8F1B0E fmul dword ptr [eax+0Ch] 0F8F1B11 faddp st(1),st 0F8F1B13 fild dword ptr [ecx+8] 0F8F1B16 fmul dword ptr [eax+10h] 0F8F1B19 faddp st(1),st 0F8F1B1B fild dword ptr [ecx+0Ch] 0F8F1B1E fmul dword ptr [eax+14h] 0F8F1B21 faddp st(1),st 0F8F1B23 fild dword ptr [ecx+10h] 0F8F1B26 fmul dword ptr [eax+18h] 0F8F1B29 faddp st(1),st 0F8F1B2B fild dword ptr [ecx+14h] 0F8F1B2E add ecx,20h 0F8F1B31 fmul dword ptr [eax+1Ch] 0F8F1B34 add eax,20h 0F8F1B37 faddp st(1),st 0F8F1B39 dec edi 0F8F1B3A jne IM_WeightedMedianBlur_03+1B5h (0F8F1AF5h) }

　　這里是明顯的普通的FPU代碼，多說一句，針對這個循環，系統也進行了多路并行優化。

　 ? 為了比較方便，我們把禁用系統優化后的時間和未禁用是做一個整體的對比：

算法名稱	執行時間
算法名稱	禁用編譯器優化	啟用編譯器優化
BruteForce	26875ms	27025ms
Joint Histgram	123432ms	108254ms
Joint Hist CacheFriend	55214ms	17325ms
Joint + MT	1075ms	420ms
Joint + NT	1286ms	1200ms
Joint + MT + NT	422ms	430ms

? ? ? 當禁用編譯器優化后，可以明顯的看到Joint + MT + NT的速度優勢比較大，和論文里給出的數據也基本相當了。

? ? ??但是我們還是稍作分析，為什么同樣是開啟編譯器優化，Joint + MT的速度能從1075ms降低到420ms，而Joint + MT + NT確基本沒有什么變化呢，這就要從代碼本身說起。

? ? ? 我們注意到，在Joint + MT版本中，BalanceWeight和CurWeight等元素的計算都是通過一個簡單的for循環進行的，計算過程中循環的次數是固定的，每次計算內部的循環變量取值也是按照內存順序來的，這種代碼非常適合編譯器使用SIMD指令優化，他會自動編譯一系列帶P（Packet）字母的SIMD指令（例如mulps）進行單周期四指令的快速執行，相當于提高了4倍的通行能力，而那些計算在整個算法里占用的時間比例有比較大，這樣對整個算法的提速表現貢獻是很大的。

? ? ? 而在有了Necklace Table參與的版本中，由于BalanceWeight和CurWeight的更新使用do while循環，循環的次數是未知的，循環里的指針指向的位置也是變動的，因此，即使使用了SIMD指令，他也只能使用其中帶S（Single）字母的SIMD指令（例如mulss），這種指令一次性也就是執行一條計算，相比普通的FPU指令提速非常有限甚至更慢，因此，優不優化速度基本沒啥區別。另外一個重要的問題在論文中其實沒有提及，那就是隨著半徑的增加，Joint Histgram中得非0元素會相對的變得越來越少（但整體比例還是很大的），但是在BCB中，只要某個固定Feature對應的LevelF個直方圖元素中有一個不為0，那么他就會不為0，這個情況在大半徑時發生的概率非常高，此時的更新Necklace Table的時間和后續減少計算的時間來說可能會本末倒置，反而會引起計算時間的增加。

　　基于這樣一個分析，隱含著這樣一個事實，當半徑比較小時，由于計算過程中非零值的存在，Joint + MT + NT應該效果會更改，而隨著半徑的增加，非零值減小，NT帶來的收益越來越小，甚至抵消了，我們實測了下面一組數據。

算法名稱	不同半徑時的執行時間（ms）
算法名稱	1	3	5	8	10	15	20	40
Joint + MT	386	404	396	416	436	500	540	744
Joint + MT + NT	153	316	306	412	452	534	654	1091

? ? ? 也就是說，在容許進行SIMD優化的情況下，當半徑大于10時，建議使用Joint + MT來獲得更高的效率，半徑小于10時，可通過Joint + MT + NT來提供更好的速度。

? ? ? 從代碼的簡練或者內存占用方面來說，毫無疑問Joint + MT更簡單，也更加節省內存，如果在現在的PC上使用該算法，我更喜歡直接使用Joint + MT算法。

? ? ? 這樣并不是說Necklace Table不好，我反到覺得這個數據結構也是由很高的利用價值，也許可以利用到我關心的其他一些算法上，會有這比較好的效果。

　　另外小聲的說一下，似乎這里的最終優化的時間和Brute force的時間比并沒有達到100:1。

? ? ? 五、后續關于Joint + MT進一步優化的幾個嘗試

　　既然選中Joint + MT，我們再仔細的構思下他還沒有進一步優化的余地呢，第一想到的就是，我自行內嵌SIMD指令，代碼中有好幾個for循環使用SIMD指令應該很容易處理，但是，經過多次改寫，發現這種非常簡便的for循環，我們自己內嵌的SIMD指令很難超越編譯器編譯后的速度，畢竟寫編譯器的那些專家的優化水平，不是我等能夠比擬的。第一步方向選擇放棄。

? ? ? 那么如果考慮定點話呢，一般兩個像素之間的權重值是個介于0和1之間的數據，如果我們把它放大一定倍數，轉換為整形，那么整個計算過程就是整形的處理，而且現在整形也可以直接使用SSE處理，同樣是一次性處理4個32位整形，同浮點相比，少了幾次數據類型的轉換，經過測試，這樣處理后速度基本沒有什么大的差異，這個方法也可以放棄。

? ? ?第三個想法是直方圖的更新，有一種常用的直方圖更新方法是特例化處理圖像整體最左上角的點，然后在水平方向移動時，去除最左側的一列信息，加上最右側的信息，當移動到第一行最右側的像素點時，此時的更新方向不是直接跳到第二行首像素，而是從第二行尾像素向第二行手像素進行處理，這時我們可以充分利用第一行的最右側像素的直方圖數據，只要減去最上部一行的直方圖信息，然后加上最下部一行的直方圖的信息就可以了，在逆向移動時，直方圖的更新則和第一行的更新相反，加上左側的信息，然后減去右側信息，當處理到第二行首地址像素后，我們又跳到第三行首地址，然后進行類似第一行的處理，這種處理方式能夠減少對每行首像素進行全部直方圖更新的計算量，在半徑較大時有一定的加速作用，我們一般稱之為蛇形算法。實驗了一下，對算法的速度提升非常有限，而且會使得代碼稍顯繁瑣。也需要放棄。

? ? ?那么目前我想到的唯一的有可能對速度還有提升的就是定點化時不用32位的數據，適當的考慮數據的范圍，如果能保證定點后的數據能在16位的有效范圍，那么還是有可能進一步提高點速度的，畢竟這個時候可以使用SSE單指令一次性進行8個整數的加減乘法了，這個有待于進一步去測試。

　　六、特例優化

　　在有些情況下甚至很多情況下，我們使用的Feature是其自身，這種情況下因為數據的特殊性，我們可以做一些特殊處理，使得算法的速度更快。

　　當Feature等于Input本身時，我們注意到，聯合直方圖中只有45度的對角線中元素有值，其他部位都為0，因此，我們可以考慮聯合直方圖在形式上退化為一維直方圖，這個時候一個簡單的代碼如下所示：

int IM_WeightedMedianBlur_Special(unsigned char *Input, float *Weight, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1) && (Channel != 3) && (Channel != 4)) return IM_STATUS_NOTSUPPORTED;const int Level = 256;int *Histgram = (int *)malloc(Level * sizeof(int));if (Histgram == NULL) return IM_STATUS_OUTOFMEMORY;for (int Y = 0; Y < Height; Y++){unsigned char *LinePS = Input + Y * Stride;unsigned char *LinePD = Output + Y * Stride;memset(Histgram, 0, Level * sizeof(int)); // 全部賦值為0for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Index = J * Stride;for (int I = IM_Max(0 - Radius, 0); I <= IM_Min(0 + Radius, Width - 1); I++){Histgram[Input[Index + I]]++;}}for (int X = 0; X < Width; X++){int IndexF = LinePS[X] * Level;float SumW = 0, HalfSumW = 0;;for (int I = 0; I < Level; I++){SumW += Histgram[I] * Weight[IndexF + I];}HalfSumW = SumW / 2;SumW = 0;for (int I = 0; I < Level; I++){SumW += Histgram[I] * Weight[IndexF + I];if (SumW >= HalfSumW){LinePD[X] = I;break;}}if ((X - Radius) >= 0){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){Histgram[Input[J * Stride + X - Radius]]--;}}if ((X + Radius + 1) <= Width - 1){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){Histgram[Input[J * Stride + X + Radius + 1]]++;}}}}free(Histgram);return IM_STATUS_OK; }

　　同樣是21*21的窗口，one - metalpixel的灰度圖像計算用367ms左右，比上述都要快。

　　同樣的道理，我們也可以使用BCB技術來優化，但是此時的BCB來的更簡單。

int IM_WeightedMedianBlur_Special_BCB(unsigned char *Input, float *Weight, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1) && (Channel != 3)) return IM_STATUS_NOTSUPPORTED;int Status = IM_STATUS_OK;const int Level = 256; int *Histgram = (int *)malloc(Level * sizeof(int));int *BCB = (int *)malloc(Level * sizeof(int));if ((Histgram == NULL) || (BCB == NULL)){Status = IM_STATUS_OK;goto FreeMemory;}for (int Y = 0; Y < Height; Y++){unsigned char *LinePS = Input + Y * Stride;unsigned char *LinePD = Output + Y * Stride;memset(Histgram, 0, Level * sizeof(int)); // 全部賦值為0memset(BCB, 0, Level * sizeof(int));int CutPoint = -1;for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Index = J * Stride;for (int I = IM_Max(0 - Radius, 0); I <= IM_Min(0 + Radius, Width - 1); I++){int Value = Input[J * Stride + I];Histgram[Value]++; // 計算每行第一個點的二維直方圖，直方圖的水平方向為Feature坐標，垂直方向為Value坐標 BCB[Value]--; // 此時的CutPoint初始化為-1，所以+方向的數據為0，所有的都在-方向 }}for (int X = 0; X < Width; X++){float BalanceWeight = 0;int IndexF = LinePS[X] * Level; // 中心點P的Value所對應的那一行Feature權重起始索引for (int I = 0; I < Level; I++) // BCB[I]中保存的是以CutPoint為分界線，Feature為I時，分界線左側的所有Value[0-CutPoint]值的數量和分界線右側所有的Value（CutPoint， LevelV - 1]值數量的差異{BalanceWeight += BCB[I] * Weight[IndexF + I]; // 因為Feature為固定值時，如果中心點固定，那么不管與Feature對應的Value值時多少，Weight就是定值了。}if (BalanceWeight < 0) // 第一個點的BalanceWeight必然小于0{for (; BalanceWeight < 0 && CutPoint != Level - 1; CutPoint++){int Index = CutPoint + 1; // 新的直方圖的位置BCB[Index] += Histgram[Index] * 2; // 數量是同樣的道理BalanceWeight += 2 * Histgram[Index] * Weight[IndexF + Index];}}else if (BalanceWeight > 0) // 如果平衡值大于0，則向左移動中間值{for (; BalanceWeight > 0 && CutPoint != 0; CutPoint--){BCB[CutPoint] -= Histgram[CutPoint] * 2;BalanceWeight -= 2 * Histgram[CutPoint] * Weight[IndexF + CutPoint];;}}LinePD[X] = CutPoint;if ((X - Radius) >= 0){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++) // 即將移出的那一列數據{int Value = Input[J * Stride + X - Radius];Histgram[Value]--;if (Value <= CutPoint) // 如果移出的那個值小于當前的中值BCB[Value]--;elseBCB[Value]++;}}if ((X + Radius + 1) <= Width - 1){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X + Radius + 1];Histgram[Value]++;if (Value <= CutPoint) // 如果移出的那個值小于當前的中值BCB[Value]++;elseBCB[Value]--;}}}} FreeMemory:if (Histgram != NULL) free(Histgram);if (BCB != NULL) free(BCB);return Status; }

　　同樣是21*21的窗口，one - metalpixel的灰度圖像計算用242ms左右。

　　如果我們進一步退化，將其退化為普通的中值濾波，即所有Weight都相同，則刪減不需要的相關代碼后，可以有如下過程：

int IM_MedianBlur(unsigned char *Input, unsigned char *Output, int Width, int Height, int Stride, int Radius) {int Channel = Stride / Width;if ((Input == NULL) || (Output == NULL)) return IM_STATUS_NULLREFRENCE;if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER;if ((Channel != 1) && (Channel != 3)) return IM_STATUS_NOTSUPPORTED;int Status = IM_STATUS_OK;const int Level = 256;int *Histgram = (int *)malloc(Level * sizeof(int));if ((Histgram == NULL)){Status = IM_STATUS_OK;goto FreeMemory;}for (int Y = 0; Y < Height; Y++){unsigned char *LinePS = Input + Y * Stride;unsigned char *LinePD = Output + Y * Stride;memset(Histgram, 0, Level * sizeof(int)); // 全部賦值為0int CutPoint = -1;int Balance = 0;for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Index = J * Stride;for (int I = IM_Max(0 - Radius, 0); I <= IM_Min(0 + Radius, Width - 1); I++){int Value = Input[J * Stride + I];Histgram[Value]++; // 計算每行第一個點的二維直方圖，直方圖的水平方向為Feature坐標，垂直方向為Value坐標 Balance--;}}for (int X = 0; X < Width; X++){ if (Balance < 0) // 第一個點的Balance必然小于0{for (; Balance < 0 && CutPoint != Level - 1; CutPoint++){ Balance += 2 * Histgram[CutPoint + 1];}}else if (Balance > 0) // 如果平衡值大于0，則向左移動中間值{for (; Balance > 0 && CutPoint != 0; CutPoint--){Balance -= 2 * Histgram[CutPoint];}}LinePD[X] = CutPoint;if ((X - Radius) >= 0){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++) // 即將移出的那一列數據{int Value = Input[J * Stride + X - Radius];Histgram[Value]--;if (Value <= CutPoint) // 如果移出的那個值小于當前的中值Balance--;elseBalance++;}}if ((X + Radius + 1) <= Width - 1){for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++){int Value = Input[J * Stride + X + Radius + 1];Histgram[Value]++;if (Value <= CutPoint) // 如果移出的那個值小于當前的中值Balance++;elseBalance--;}}}} FreeMemory:if (Histgram != NULL) free(Histgram);return Status; }

? ? ?同樣是21*21的窗口，one - metalpixel的灰度圖像計算用140ms左右。

? ? ?有興趣的朋友還可以試下對上述中值模糊的代碼在加上Necklace table優化，看看能得到什么樣的結果。

? ? ?在論文的最后，講述了加權中值模糊的多個應用場景，比如在光流、立體匹配、JPG瑕疵修復、藝術特效等等方面，我測試下幾個我能做的測試，確實有不錯的效果，比如下面的JPG瑕疵修復。對簡單的圖處理后確實蠻好的，如果在結合我之前研究的MLAA去鋸齒算法，恢復后的圖像質量就更高了，如下所示：

? ? ??? ? ? ? ??? ? ? ? ? ?

? ?? ? ?? ? ?

? ? ? ? ?? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ?

　　　　　　　　　　　　原圖　　　　　　　　　　　　　　　　　　　　　　　　　　 ? ? ?　加權中值模糊（特征圖為原圖）　　　　　　　　　　　　　　　　　 MLAA后續處理后（邊緣更平滑）

? ? ? ? 另外，WMF的保邊特性感覺比其他的如導向濾波、雙邊濾波等等都要強烈的多，比如下圖： ?

? ?? ? ??

　　花朵的邊緣，下面的文字等等處理都還特別清晰，不像其他的保邊濾波器總有點模糊，這個特性也許用到一些增強上也會有很不錯的效果。

　　按照上述文章的思路，我整理和編制一個簡易的測試程序，用來論證論文和我博文中得一些數據，使用的VS2013編譯的，用C++做的DLL，C#做的UI測試界面，不依賴于任何其他第三方庫，目前只做了灰度圖的方案，因為彩色的話也基本就是三個通道獨立寫，可以通過拆分然后調用灰度的來實現。我也測試了下作者分享的VS工程，應該比我提供的代碼速度稍微慢一點。

　　源代碼下載地址：https://files.cnblogs.com/files/Imageshop/WeightedMedianBlur.rar

總結

以上是生活随笔為你收集整理的加权中值滤波的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

中值