CVPR論文《100+ Times FasterWeighted Median Filter (WMF)》的實現和解析（附原始碼）。 任意半徑中值濾波（擴充套件至百分比濾波器）O(1)時間複雜度演算法的原理、實現及效果 任意半徑中值濾波（擴充套件至百分比濾波器）O(1)時間複雜度演算法的原理、實現及
阿新 • • 發佈：2018-11-13
四年前第一次看到《100+ Times FasterWeighted Median Filter (WMF)》一文時，因為他附帶了原始碼，而且還是CVPR論文，因此，當時也對程式碼進行了一定的整理和解讀，但是當時覺得這個演算法雖然對原始速度有不少的提高，但是還是比較慢。因此，沒有怎麼在意，這幾天有幾位朋友又提到這篇文章，於是把當時的程式碼和論文又仔細的研讀了一番，對論文的思想和其中的實現也有了一些新的新的，再次做個總結和分享。
這篇文章的官網地址是：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;
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的效能分析器，可以大概獲得如下的結果：
可見核心程式碼基本都用於排序了，使用更快的排序有助於進一步提高速度。
針對這個情況，論文的作者從多方面提出了改進措施，主要有三個方面，我們簡單的重複下。
一、聯合直方圖(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)); // 全部賦值為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)。提速也是相當的可觀，而且演算法速度和半徑不是特別敏感，畢竟更新直方圖的計算量在這裡佔的比例其實已經不多了。
三、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倍，由於這個演算法和影象內容是由一定關係的，因此，和論文提供的資料直接比較的意義不大。
四、最終的結合體
很自然的，我們想到要把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]值數量的差異 { BalanceWei