CVPR論文《100+ Times Faster Weighted Median Filter (WMF)》的實現和解析（附源代碼）。
四年前第一次看到《100+ Times FasterWeighted Median Filter (WMF)》一文時，由於他附帶了源代碼，並且仍是CVPR論文，所以，當時也對代碼進行了必定的整理和解讀，可是當時以爲這個算法雖然對原始速度有很多的提升，可是仍是比較慢。所以，沒有怎麼在乎，這幾天有幾位朋友又提到這篇文章，因而把當時的代碼和論文又仔細的研讀了一番，對論文的思想和其中的實現也有了一些新的新的，再次作個總結和分享。html
這篇文章的官網地址是：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。數據結構
因爲做者只給出了最後的優化實現代碼，而論文中還提出了各類中間過程的時間，所以本文以實現和驗證論文中有關說法爲主，涉及到的理論知識比較膚淺，通常是一筆而過。ide
根據論文中得說法，所謂的加權中值濾波，也是一種非線性的圖像平滑技術，他取一個局部窗口內全部像素的加權中值來代替局部窗口的中心點的值。用較爲數學的方法表示以下：函數
在圖像I中的像素p，咱們考慮以p爲中心，半徑爲R的局部窗口
，不一樣於普通的中值模糊，對於屬於內每個像素q，都有一個基於對應的特徵圖像的類似度的權重係數wpq，以下式所示：佈局
f(p)和f(q)是像素p和q在對應的特徵圖中得特徵值。g是一個權重函數，最經常使用的即爲高斯函數，反應了像素p和q的類似程度。post
咱們用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; else
return 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的性能分析器，能夠大概得到以下的結果：
可見核心代碼基本都用於排序了，使用更快的排序有助於進一步提升速度。
針對這個狀況，論文的做者從多方面提出了改進措施，主要有三個方面，咱們簡單的重複下。
1、聯合直方圖(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），所以實際上，中值查找這一塊的耗時佔了絕對的比例。
2、快速中值追蹤
尋找中值的過程實際上能夠當作一個追求平衡的過程，假定當前搜索到的位置是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)); // 所有賦值爲0 memset(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]--; else BCB[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]++; else BCB[Feature]--; } } } } free(Histgram); free(BCB); }
代碼也很簡潔，主要是增長了一個BCB列表的維護，編譯後測試，一樣是21*21的窗口，one - metalpixel的灰度圖像計算用420ms左右，比Brute-force版本的27s大約快了64倍，這個和論文的時間比例基本差很少(（156.9+0.4）/(2.2+0.5)=58)。提速也是至關的可觀，並且算法速度和半徑不是特別敏感，畢竟更新直方圖的計算量在這裏佔的比例其實已經很少了。
3、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 table int *BackWordH = (int *)malloc(LevelF * LevelV * sizeof(int)); // forward link for necklace table float *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，須要把第一個元素置爲0 BackWordH[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倍，因爲這個算法和圖像內容是由必定關係的，所以，和論文提供的數據直接比較的意義不大。
4、最終的結合體
很天然的，咱們想到要把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 table int *BackWordH = (int *)malloc(LevelF * LevelV * sizeof(int)); // forward link for necklace table int *ForwardBCB = (int *)malloc(LevelF * sizeof(int)); // forward link for necklace table int *BackWordBCB = (int *)malloc(LevelF * sizeof(int)); // forward link for necklace table if ((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)); // 所有賦值爲0 memset(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的算法的下面一句代碼看看：
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就是定值了。 }
反編譯結果爲：
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就是定值了。 0FAF1B25 movdqu xmm0,xmmword ptr [ecx] 0FAF1B29 add ecx,10h 0FAF1B2C cvtdq2ps xmm1,xmm0 0FAF1B2F movups xmm0,xmmword ptr [eax] 0FAF1B32 add eax,10h 0FAF1B35 mulps xmm1,xmm0 0FAF1B38 addps xmm2,xmm1 0FAF1B3B dec edx 0FAF1B3C jne IM_WeightedMedianBlur_03+1B5h (0FAF1B25h) }
赤裸裸的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。
5、後續關於Joint + MT進一步優化的幾個嘗試
既然選中Joint + MT，咱們再仔細的構思下他尚未進一步優化的餘地呢，第一想到的就是，我自行內嵌SIMD指令，代碼中有好幾個for循環使用SIMD指令應該很容易處理，可是，通過屢次改寫，發現這種很是簡便的for循環，咱們本身內嵌的SIMD指令很難超越編譯器編譯後的速度，畢竟寫編譯器的那些專家的優化水平，不是我等可以比擬的。第一步方向選擇放棄。
那麼若是考慮定點話呢，通常兩個像素之間的權重值是個介於0和1之間的數據，若是咱們把它放大必定倍數，轉換爲整形，那麼整個計算過程就是整形的處理，並且如今整形也能夠直接使用SSE處理，一樣是一次性處理4個32位整形，同浮點相比，少了幾回數據類型的轉換，通過測試，這樣處理後速度基本沒有什麼大的差別，這個方法也能夠放棄。
第三個想法是直方圖的更新，有一種經常使用的直方圖更新方法是特例化處理圖像總體最左上角的點，而後在水平方向移動時，去除最左側的一列信息，加上最右側的信息，當移動到第一行最右側的像素點時，此時的更新方向不是直接跳到第二行首像素，而是從第二行尾像素向第二行手像素進行處理，這時咱們能夠充分利用第一行的最右側像素的直方圖數據，只要減去最上部一行的直方圖信息，而後加上最下部一行的直方圖的信息就能夠了，在逆向移動時，直方圖的更新則和第一行的更新相反，加上左側的信息，而後減去右側信息，當處理到第二行首地址像素後，咱們又跳到第三行首地址，而後進行相似第一行的處理，這種處理方式可以減小對每行首像素進行所有直方圖更新的計算量，在半徑較大時有必定的加速做用，咱們通常稱之爲蛇形算法。實驗了一下，對算法的速度提高很是有限，並且會使得代碼稍顯繁瑣。也須要放棄。
那麼目前我想到的惟一的有可能對速度還有提高的就是定點化時不用32位的數據，適當的考慮數據的範圍，若是能保證定點後的數據能在16位的有效範圍，那麼仍是有可能進一步提升點速度的，畢竟這個時候可使用SSE單指令一次性進行8個整數的加減乘法了，這個有待於進一步去測試。
6、特例優化
在有些狀況下甚至不少狀況下，咱們使用的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)); // 所有賦值爲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++) { 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)); // 所有賦值爲0 memset(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]--; else BCB[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]++; else BCB[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)); // 所有賦值爲0 int 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--; else Balance++; } } 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++; else Balance--; } } } } 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
後記：
寫完文章後，對於Joint + MT算法總以爲應該還有能夠繼續改進的地方，這幾日也還在琢磨這事，有如下幾個收穫能夠進一步提升速度。
第一：正如前文所描述的是否能夠考慮直方圖數據用16位來表示呢，包括BCB都用16位。咱們來簡單分析下。
咱們看到Joint Histgram中一共有LevelV * LevelF個元素，局部窗口內有n=(2r+1)*(2r+1)個元素，那麼在最極端的狀況下這n個元素值都相同，且其對應的n個feature值也相同，這樣在histgram中元素的最大值就是n，只要這個n小於short能表示的最大的正數，則用short就能夠徹底表達完整的直方圖信息，此時對應的r值約爲90，徹底能知足實際的需求了，並且這種極端狀況基本不會發生。
一樣對於BCB，也不太可能在一個局部出現其值超出short能表示的正負範圍的。
那麼對於權重值，咱們也能夠把他們定點化，通常權重都會在[0,1]範圍內的一個數，即便不是，咱們也能夠把他們歸一化，而後好比放大16384倍，使用一個short類型數據來保存。
這樣作的好處就是，咱們可使用simd中關於16位的一些高級計算指令了，好比下面這段代碼：
int BalanceWeight = 0; for (int I = 0; I < LevelF; I++) { BalanceWeight += BCB[I] * W[IndexF + I]; }
則能夠優化爲：
__m128i BW1 = _mm_setzero_si128(); __m128i BW2 = _mm_setzero_si128(); for (int I = 0; I < Block * BlockSize; I += BlockSize) { BW1 = _mm_add_epi32(BW1, _mm_madd_epi16(_mm_load_si128((__m128i *)(BCB + I)), _mm_load_si128((__m128i *)(W + IndexF + I)))); BW2 = _mm_add_epi32(BW2, _mm_madd_epi16(_mm_load_si128((__m128i *)(BCB + I + 8)), _mm_load_si128((__m128i *)(W + IndexF + I + 8)))); } int BalanceWeight = _mm_hsum_epi32(_mm_add_epi32(BW1, BW2));
其中int BlockSize = 16, Block = LevelF / BlockSize;
_mm_madd_epi16能夠一次性的進行8個16位數據的乘法和加法計算，效率及其高效。
後面的BalanceWeight 中得校訂代碼的for循環也可使用相似方法優化。
這樣作還有個好處就是佔用的內存小了，並且Y循環裏的memset工做量也會少一半。
第2、咱們在更新BCB的時候一段這樣的小代碼：
if (Value <= CutPoint) BCB[Feature]--; else BCB[Feature]++;
別看很短小，因爲他出如今直方圖的更新裏，所以執行的頻率很高。咱們看下他的反彙編：
if (Value <= CutPoint) 0F271E3D mov edi,dword ptr [esp+28h] 0F271E41 cmp ecx,dword ptr [esp+0Ch] 0F271E45 jg IM_WeightedMedianBlur_Joint_MT+495h (0F271E55h) BCB[Feature]--; 0F271E47 mov ecx,dword ptr [esp+38h] 0F271E4B dec word ptr [edi+edx*2] Index += Stride; 0F271E4F add esi,dword ptr [Stride] 0F271E52 inc ecx 0F271E53 jmp IM_WeightedMedianBlur_Joint_MT+44Dh (0F271E0Dh) for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++) // 即將移出的那一列數據 0F271E55 mov ecx,dword ptr [esp+38h] else BCB[Feature]++; 0F271E59 inc word ptr [edi+edx*2] Index += Stride; 0F271E5D add esi,dword ptr [Stride] 0F271E60 inc ecx 0F271E61 jmp IM_WeightedMedianBlur_Joint_MT+44Dh (0F271E0Dh) Index += Stride;
很明顯，裏面有jmp跳轉指令，之前我不以爲這個有什麼速度影響，可是咱們嘗試着用一些其餘代碼技巧替代這段代碼後，速度卻有了質的提高。
BCB[Value] += -((Value <= CutPoint) << 1) + 1;
反彙編看下：
BCB[Value] += -((Value <= CutPoint) << 1) + 1; 0F93126B mov eax,1 0F931270 mov ecx,dword ptr [ebp-20h] 0F931273 dec word ptr [ecx+edx*2] 0F931277 xor ecx,ecx 0F931279 cmp edx,esi 0F93127B setle cl 0F93127E add cx,cx 0F931281 sub ax,cx 0F931284 mov ecx,dword ptr [ebp-40h] 0F931287 add word ptr [edi+edx*2],ax 0F93128B mov eax,dword ptr [ebp-14h] 0F93128E inc eax 0F93128F mov dword ptr [ebp-14h],eax
裏面沒有了jmp跳轉了。
第三：程序裏的IM_Max有點多，能夠提取到X的循環外面。
第四個嘗試是，咱們在更新直方圖時是按列更新的，這種狀況的Cache Miss至關嚴重，一種改進的方式是，咱們備份一個原圖和特徵圖的轉置圖，這個時候更新直方圖時就是按照行方向讀取數據了，此時會多一個轉置的操做，可是轉置已經可使用SSE代碼進行高度的優化，在這個算法裏這個的耗時幾乎能夠忽略不計。實測在半徑小於10時，對總體速度無啥影響，可是隨着半徑增大，這個的效果愈來愈明顯了。
通過上述幾個步驟的優化和處理，一樣是21*21的窗口，one - metalpixel的灰度圖像計算用時205ms左右，若是是彩色圖像耗時大概在630ms左右，這個比做者提供的代碼的執行速度快了大概3倍。
在Input和Feature相同的狀況下，也能夠作一樣的優化，此時，一樣是21*21的窗口，one - metalpixel的灰度圖像計算用時125ms左右，若是是彩色圖像耗時大概在340ms左右。
一個示例Demo可從：SSE_Optimization_Demo.rar處下載。
若是以爲本文對你有幫助，還請點個贊或者給我來杯小Coffee吧。
- 1. Faster RCNN算法解析（附源代碼，能夠直接運行）
- 2. CVPR 2020論文開源項目一頁看盡，附代碼論文
- 3. CVPR 2020 論文開源項目一頁看盡，附代碼論文
- 4. CVPR 2014 論文
- 5. Faster R-CNN論文及源碼解讀
- 6. CVPR/AAAI/ECCV頂會論文/代碼
- 7. faster rcnn源碼解析
- 8. java解析WMF文件
- 9. 詳細的Faster R-CNN源碼解析之RPN源碼解析
- 10. CVPR 2020代碼開源的論文最全合集
- 更多相關文章...
- • TCP報文格式解析 - TCP/IP教程
- • Markdown 代碼 - Markdown 教程
- • Scala 中文亂碼解決
- • ☆基於Java Instrument的Agent實現
-
每一个你不满意的现在，都有一个你没有努力的曾经。