當(dāng)前位置：首頁(yè) > 编程资源 > 编程问答 >内容正文

编程问答

加权中值滤波

發(fā)布時(shí)間：2025/3/17 编程问答 24 豆豆

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

應(yīng)作者要求，轉(zhuǎn)載自：https://www.cnblogs.com/Imageshop/p/9934670.html

《100+ Times FasterWeighted Median Filter (WMF)》

　　這篇文章的官網(wǎng)地址是：http://www.cse.cuhk.edu.hk/~leojia/projects/fastwmedian/，其中主要作者Jiaya Jia教授的官網(wǎng)地址是：http://jiaya.me/，根據(jù)Jiaya Jia的說(shuō)法，這個(gè)算法很快將被OpenCv所收錄，到時(shí)候OpenCv的大神應(yīng)該對(duì)他還有所改進(jìn)吧。

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

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

　　根據(jù)論文中得說(shuō)法，所謂的加權(quán)中值濾波，也是一種非線性的圖像平滑技術(shù)，他取一個(gè)局部窗口內(nèi)所有像素的加權(quán)中值來(lái)代替局部窗口的中心點(diǎn)的值。用較為數(shù)學(xué)的方法表示如下：

　　在圖像I中的像素p，我們考慮以p為中心，半徑為R的局部窗口，不同于普通的中值模糊，對(duì)于屬于內(nèi)每一個(gè)像素q，都有一個(gè)基于對(duì)應(yīng)的特征圖像的相似度的權(quán)重系數(shù)wpq，如下式所示：

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

　　f(p)和f(q)是像素p和q在對(duì)應(yīng)的特征圖中得特征值。g是一個(gè)權(quán)重函數(shù)，最常用的即為高斯函數(shù)，反應(yīng)了像素p和q的相似程度。

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

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

　　很明顯，上面的過(guò)程要比標(biāo)準(zhǔn)的中值模糊復(fù)雜一些，在處理時(shí)多了特征圖和權(quán)重函數(shù)項(xiàng)，而標(biāo)準(zhǔn)的中值模糊我們可以認(rèn)為是加權(quán)中值模糊的特例，即所有局部窗口的權(quán)重都為1或者說(shuō)相等。

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

　　按照這個(gè)定義，我們給出一段針對(duì)灰度數(shù)據(jù)的Brute-force處理代碼：

int __cdecl ComparisonFunction(const void *X, const void *Y) // 一定要用__cdecl這個(gè)標(biāo)識(shí)符 {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; }// 加權(quán)中值模糊，直接按照算法的定義實(shí)現(xiàn)。 // Input - 輸入圖像，灰度圖，LevelV = 256級(jí) // FeatureMap - 特征圖像，灰度圖，LevelF = 256級(jí) // Weight - 特征的權(quán)重矩陣，大小是LevelF * LevelF // Output - 輸出圖像，不能和Input為同一個(gè)數(shù)據(jù)。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 可能出現(xiàn)的不同數(shù)量const int LevelF = 256; // Feature 可能出現(xiàn)的不同數(shù)量Value_Weight *VW = (Value_Weight *)malloc((2 * Radius + 1) * (2 * Radius + 1) * sizeof(Value_Weight)); // 值和特征序列對(duì)內(nèi)存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]; // 對(duì)應(yīng)的權(quán)重VW[PixelAmount].Value = Value;VW[PixelAmount].Weight = CurWeight; // 保存數(shù)據(jù)SumW += CurWeight; // 計(jì)算累加數(shù)據(jù)PixelAmount++; // 有效的數(shù)據(jù)量 }}float HalfSumW = SumW * 0.5f; // 一半的權(quán)重SumW = 0;qsort(VW, PixelAmount, sizeof VW[0], &ComparisonFunction); // 調(diào)用系統(tǒng)的qsort按照Value的值從小到大排序，注意qsort的結(jié)果仍然保存在第一個(gè)參數(shù)中for (int I = 0; I < PixelAmount; I++) // 計(jì)算中值{SumW += VW[I].Weight;if (SumW >= HalfSumW){LinePD[X] = VW[I].Value;break;}}}}free(VW);return IM_STATUS_OK; }

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

　　我們?cè)谝慌_(tái) I5，3.3GHZ的機(jī)器上進(jìn)行了測(cè)試，上述代碼處理一副1000*1000像素的灰度圖，半徑為10（窗口大小21*21）時(shí)，處理時(shí)間約為27s，論文里給的Cpu和我的差不多，給出的處理one - metalpixel的RGB圖用時(shí)90.7s，考慮到RGB的通道的數(shù)據(jù)量以及一些其他的處理，應(yīng)該說(shuō)論文如實(shí)匯報(bào)了測(cè)試數(shù)據(jù)。

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

? ? ? ?

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

　　針對(duì)這個(gè)情況，論文的作者從多方面提出了改進(jìn)措施，主要有三個(gè)方面，我們簡(jiǎn)單的重復(fù)下。

　　一、聯(lián)合直方圖(Joint Histgram）

　　直方圖優(yōu)化在很多算法中都有應(yīng)用，比如標(biāo)準(zhǔn)的中值濾波，現(xiàn)在看到的最快的實(shí)現(xiàn)方式還是基于直方圖的，詳見(jiàn)：任意半徑中值濾波（擴(kuò)展至百分比濾波器）O(1)時(shí)間復(fù)雜度算法的原理、實(shí)現(xiàn)及效果，但是在加權(quán)中值濾波中，傳統(tǒng)的一維直方圖已經(jīng)無(wú)法應(yīng)用，因?yàn)檫@個(gè)算法不僅涉及到原圖的像素值，還和另外一幅特征圖有關(guān)，因此，文中提出了聯(lián)合直方圖，也是一種二維直方圖。

　　如果圖像中的像素最多有LevelV個(gè)不同值，其對(duì)應(yīng)的特征最多有LevelF個(gè)不同的值，那么我們定義一個(gè)寬和高分別為L(zhǎng)evelV * LevelF大小的直方圖。對(duì)于某一個(gè)窗口，統(tǒng)計(jì)其內(nèi)部的(2r+1)*(2r+1)個(gè)像素和特征對(duì)的直方圖數(shù)據(jù)，即如果某個(gè)點(diǎn)的像素值為V，對(duì)應(yīng)的特征值為F，則相應(yīng)位置的直方圖數(shù)據(jù)加1。

　　如果我們統(tǒng)計(jì)出這個(gè)二維的直方圖數(shù)據(jù)后，由于中心點(diǎn)的特征值是固定的，因此，對(duì)于直方圖的每一個(gè)LevelF值，權(quán)重是一定的了，我們只需計(jì)算出直方圖內(nèi)每一個(gè)Value值所對(duì)應(yīng)所有的Feature的權(quán)重后，就可方便的統(tǒng)計(jì)出中值所在的位置了。

　　那么如果每個(gè)像素點(diǎn)都進(jìn)行領(lǐng)域直方圖的計(jì)算，這個(gè)的工作量也是蠻大的，同一維直方圖的優(yōu)化思路一樣，在進(jìn)行逐像素行處理的時(shí)候，對(duì)直方圖數(shù)據(jù)可以進(jìn)行逐步的更新，去除掉移走的那一列的直方圖信息，在加入即將進(jìn)入那一列數(shù)據(jù)，而中間重疊部分則不需要調(diào)整。

　　按照論文中的Joint Histgram的布局，即行方向大小為L(zhǎng)evelV，列方向大小為L(zhǎng)evelF，編制Joint Histgram實(shí)現(xiàn)的加權(quán)中值算法代碼如下所示：

// 加權(quán)中值模糊，基于論文中圖示的內(nèi)存布局設(shè)置的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 可能出現(xiàn)的不同數(shù)量const int LevelF = 256; // Feature 可能出現(xiàn)的不同數(shù)量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]; // 統(tǒng)計(jì)二維直方圖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++) // 計(jì)算每個(gè)Value列針對(duì)的不同的Feature的權(quán)重的累計(jì)值{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) // 計(jì)算中值{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; }

　　編譯后測(cè)試，同樣是21*21的窗口，one - metalpixel的灰度圖像計(jì)算用時(shí)多達(dá)108s，比直接實(shí)現(xiàn)慢很多了。

　　分析原因，核心就是在中值的查找上，由于我們采用的內(nèi)存布局方式，導(dǎo)致計(jì)算每個(gè)Value對(duì)應(yīng)的權(quán)重累加存在的大量的Cache miss現(xiàn)象，即下面這條語(yǔ)句：

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

　　我們換種Joint Histgram的布局，即行方向大小為L(zhǎng)evelF，列方向大小為L(zhǎng)evelV，此時(shí)的代碼如下：

// 加權(quán)中值模糊，修改內(nèi)存布局設(shè)置的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 可能出現(xiàn)的不同數(shù)量const int LevelF = 256; // Feature 可能出現(xiàn)的不同數(shù)量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; }

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

　　和原論文比較，原論文的joint histgram時(shí)間要比直接實(shí)現(xiàn)慢（156.9s vs 90.7s)，而我這里的一個(gè)版本比brute force的快，一個(gè)比brute force的慢，因此，不清楚作者在比較時(shí)采用了何種編碼方式，但是這都不重要，因?yàn)樗麄兊膮^(qū)別都還在一個(gè)數(shù)量級(jí)上。

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

　二、快速中值追蹤

　　尋找中值的過(guò)程實(shí)際上可以看成一個(gè)追求平衡的過(guò)程，假定當(dāng)前搜索到的位置是V，位于V左側(cè)所有相關(guān)值的和是Wl，位于V右側(cè)所有相關(guān)值得和是Wr，則中值的尋找可以認(rèn)為是下式：

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

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

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

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

// 加權(quán)中值模糊, 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 可能出現(xiàn)的不同數(shù)量const int LevelF = 256; // Feature 可能出現(xiàn)的不同數(shù)量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]++; // 計(jì)算每行第一個(gè)點(diǎn)的二維直方圖，直方圖的水平方向?yàn)镕eature坐標(biāo)，垂直方向?yàn)閂alue坐標(biāo) BCB[Feature]--; // 此時(shí)的CutPoint初始化為-1，所以+方向的數(shù)據(jù)為0，所有的都在-方向 }}for (int X = 0; X < Width; X++){float BalanceWeight = 0;int IndexF = LinePF[X] * LevelF; // 中心點(diǎn)P的Value所對(duì)應(yīng)的那一行Feature權(quán)重起始索引for (int I = 0; I < LevelF; I++) // BCB[I]中保存的是以CutPoint為分界線，Feature為I時(shí)，分界線左側(cè)的所有Value[0-CutPoint]值的數(shù)量和分界線右側(cè)所有的Value（CutPoint， LevelV - 1]值數(shù)量的差異{BalanceWeight += BCB[I] * Weight[IndexF + I]; // 因?yàn)镕eature為固定值時(shí)，如果中心點(diǎn)固定，那么不管與Feature對(duì)應(yīng)的Value值時(shí)多少，Weight就是定值了。}if (BalanceWeight < 0) // 第一個(gè)點(diǎn)的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]; // 左側(cè)加右側(cè)同時(shí)減，所以是2倍BCB[I] += Histgram[IndexH + I] * 2; // 數(shù)量是同樣的道理}BalanceWeight += CurWeight;}}else if (BalanceWeight > 0) // 如果平衡值大于0，則向左移動(dòng)中間值{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++) // 即將移出的那一列數(shù)據(jù){int Value = Input[J * Stride + X - Radius];int Feature = FeatureMap[J * Stride + X - Radius];Histgram[Value * LevelF + Feature]--;if (Value <= CutPoint) // 如果移出的那個(gè)值小于當(dāng)前的中值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) // 如果移出的那個(gè)值小于當(dāng)前的中值BCB[Feature]++;elseBCB[Feature]--;}}}}free(Histgram);free(BCB); }

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

　　三、Necklace Table

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

? ??

　　比如上圖是某個(gè)圖像局部窗口的聯(lián)合直方圖和BCB值，在聯(lián)合直方圖中大部分區(qū)域都是0值對(duì)應(yīng)的黑色，在BCB中大部分情況也是0值。

? ? ? ?因此，作者構(gòu)建了一個(gè)叫做Necklace Table的數(shù)據(jù)結(jié)構(gòu)，這個(gè)數(shù)據(jù)結(jié)構(gòu)可以方便快捷的記錄下一個(gè)和上一個(gè)非0元素的位置，從而能有效的訪問(wèn)到那些真正有計(jì)算價(jià)值的部位，以及簡(jiǎn)單的刪除和增加節(jié)點(diǎn)的功能，具體的實(shí)現(xiàn)細(xì)節(jié)詳見(jiàn)論文或下面的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 可能出現(xiàn)的不同數(shù)量const int LevelF = 256; // Feature 可能出現(xiàn)的不同數(shù)量 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; // 其實(shí)每一個(gè)Feature對(duì)應(yīng)一個(gè)完整的Necklace Table，需要把第一個(gè)元素置為0BackWordH[X * LevelF] = 0;}for (int J = IM_Max(Y - Radius, 0); J <= IM_Min(Y + Radius, Height - 1); J++) // 第一個(gè)元素{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) // 直方圖數(shù)據(jù)如果還是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]; // 跳過(guò)那些非0的元素J = ForwardH[Index + J];} while (J != 0);Sum[I] = Cum; // 計(jì)算每一個(gè)Value對(duì)應(yīng)的所有Featrue的權(quán)重累計(jì)和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) // 直方圖數(shù)據(jù)如果還是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; }

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

? ? ?四、最終的結(jié)合體

　很自然的，我們想到要把Median Tracking 和 Necklace Table聯(lián)合在一起，來(lái)進(jìn)一步的提高速度，這個(gè)時(shí)候可以對(duì)Joint Histgram即BCB都使用?Necklace Table來(lái)記錄非零元素，于是產(chǎn)生了以下的結(jié)合代碼：

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) // 直方圖數(shù)據(jù)如果還是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]++; // 計(jì)算每行第一個(gè)點(diǎn)的二維直方圖，直方圖的水平方向?yàn)镕eature坐標(biāo)，垂直方向?yàn)閂alue坐標(biāo) UpdateBCB(BCB[Feature], ForwardBCB, BackWordBCB, Feature, -1); // 此時(shí)的CutPoint初始化為-1，所以+方向的數(shù)據(jù)為0，所有的都在-方向 }}for (int X = 0; X < Width; X++){float BalanceWeight = 0;int IndexF = LinePF[X] * LevelF; // 中心點(diǎn)P的Value所對(duì)應(yīng)的那一行Feature權(quán)重起始索引int I = 0;do{BalanceWeight += BCB[I] * Weight[IndexF + I]; // 按照當(dāng)前BCB數(shù)據(jù)計(jì)算平衡值，BCB記錄了相同的FMap值時(shí)按照之前的中間值左右兩側(cè)像素個(gè)數(shù)的差異值I = ForwardBCB[I];} while (I != 0);if (BalanceWeight < 0) // 第一個(gè)點(diǎn)的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]; // 左側(cè)加右側(cè)同時(shí)減，所以是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，則向左移動(dòng)中間值{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]; // 左側(cè)加右側(cè)同時(shí)減，所以是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++) // 即將移出的那一列數(shù)據(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) // 直方圖數(shù)據(jù)如果還是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; }

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

　　我們仔細(xì)檢查論文里，在Implementation Notes節(jié)里有這樣的語(yǔ)句：

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

　　第一，他也是用的C++和我一樣，第二，他是單線程，也和我一樣，第三，沒(méi)有使用任何SIMD指令，似乎我也沒(méi)有使用啊，都一樣，為什么結(jié)果比對(duì)不一致，難道是大神他們作弊，鑒于他們的成就，我立即撤回我這逆天的想法，一定是其他地方有問(wèn)題。我們?cè)囍淳幾g看看。

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

　　反編譯結(jié)果為：

　　赤裸裸的SIMD指令啊。

　　為什么呢，只是因?yàn)閂S的編譯器即使在默認(rèn)情況下的設(shè)置中，也會(huì)根據(jù)當(dāng)前編譯系統(tǒng)的情況，進(jìn)行一定的向量化優(yōu)化，加上現(xiàn)在的PC基本沒(méi)有哪一個(gè)不能使用SIMD指令的。如下圖所示，為C++默認(rèn)編譯選項(xiàng)：

? ? ? ? ? ?

　　在啟用增強(qiáng)指令集選項(xiàng)里默認(rèn)是未設(shè)置，但是未設(shè)置并不代表不使用，正如上述所言，測(cè)試編譯器會(huì)根據(jù)系統(tǒng)狀況優(yōu)化編譯。因此，雖然表面上代碼沒(méi)有使用SIMD指令，但是實(shí)際卻使用了。

　　為了公平起見(jiàn)，我們禁用系統(tǒng)的SIMD優(yōu)化，此時(shí)，可以在增強(qiáng)指令集的選項(xiàng)里選擇“無(wú)增強(qiáng)指令/arch:IA32".

? ? ? ? ?

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

for (int I = 0; I < LevelF; I++) // BCB[I]中保存的是以CutPoint為分界線，Feature為I時(shí)，分界線左側(cè)的所有Value[0-CutPoint]值的數(shù)量和分界線右側(cè)所有的Value（CutPoint， LevelV - 1]值數(shù)量的差異{BalanceWeight += BCB[I] * Weight[IndexF + I]; // 因?yàn)镕eature為固定值時(shí)，如果中心點(diǎn)固定，那么不管與Feature對(duì)應(yīng)的Value值時(shí)多少，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代碼，多說(shuō)一句，針對(duì)這個(gè)循環(huán)，系統(tǒng)也進(jìn)行了多路并行優(yōu)化。

　 ? 為了比較方便，我們把禁用系統(tǒng)優(yōu)化后的時(shí)間和未禁用是做一個(gè)整體的對(duì)比：

算法名稱(chēng)	執(zhí)行時(shí)間
算法名稱(chēng)	禁用編譯器優(yōu)化	啟用編譯器優(yōu)化
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

? ? ? 當(dāng)禁用編譯器優(yōu)化后，可以明顯的看到Joint + MT + NT的速度優(yōu)勢(shì)比較大，和論文里給出的數(shù)據(jù)也基本相當(dāng)了。

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

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

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

　　基于這樣一個(gè)分析，隱含著這樣一個(gè)事實(shí)，當(dāng)半徑比較小時(shí)，由于計(jì)算過(guò)程中非零值的存在，Joint + MT + NT應(yīng)該效果會(huì)更改，而隨著半徑的增加，非零值減小，NT帶來(lái)的收益越來(lái)越小，甚至抵消了，我們實(shí)測(cè)了下面一組數(shù)據(jù)。

算法名稱(chēng)	不同半徑時(shí)的執(zhí)行時(shí)間（ms）
算法名稱(chēng)	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

? ? ? 也就是說(shuō)，在容許進(jìn)行SIMD優(yōu)化的情況下，當(dāng)半徑大于10時(shí)，建議使用Joint + MT來(lái)獲得更高的效率，半徑小于10時(shí)，可通過(guò)Joint + MT + NT來(lái)提供更好的速度。

? ? ? 從代碼的簡(jiǎn)練或者內(nèi)存占用方面來(lái)說(shuō)，毫無(wú)疑問(wèn)Joint + MT更簡(jiǎn)單，也更加節(jié)省內(nèi)存，如果在現(xiàn)在的PC上使用該算法，我更喜歡直接使用Joint + MT算法。

? ? ? 這樣并不是說(shuō)Necklace Table不好，我反到覺(jué)得這個(gè)數(shù)據(jù)結(jié)構(gòu)也是由很高的利用價(jià)值，也許可以利用到我關(guān)心的其他一些算法上，會(huì)有這比較好的效果。

　　另外小聲的說(shuō)一下，似乎這里的最終優(yōu)化的時(shí)間和Brute force的時(shí)間比并沒(méi)有達(dá)到100:1。

? ? ? 五、后續(xù)關(guān)于Joint + MT進(jìn)一步優(yōu)化的幾個(gè)嘗試

　　既然選中Joint + MT，我們?cè)僮屑?xì)的構(gòu)思下他還沒(méi)有進(jìn)一步優(yōu)化的余地呢，第一想到的就是，我自行內(nèi)嵌SIMD指令，代碼中有好幾個(gè)for循環(huán)使用SIMD指令應(yīng)該很容易處理，但是，經(jīng)過(guò)多次改寫(xiě)，發(fā)現(xiàn)這種非常簡(jiǎn)便的for循環(huán)，我們自己內(nèi)嵌的SIMD指令很難超越編譯器編譯后的速度，畢竟寫(xiě)編譯器的那些專(zhuān)家的優(yōu)化水平，不是我等能夠比擬的。第一步方向選擇放棄。

? ? ? 那么如果考慮定點(diǎn)話(huà)呢，一般兩個(gè)像素之間的權(quán)重值是個(gè)介于0和1之間的數(shù)據(jù)，如果我們把它放大一定倍數(shù)，轉(zhuǎn)換為整形，那么整個(gè)計(jì)算過(guò)程就是整形的處理，而且現(xiàn)在整形也可以直接使用SSE處理，同樣是一次性處理4個(gè)32位整形，同浮點(diǎn)相比，少了幾次數(shù)據(jù)類(lèi)型的轉(zhuǎn)換，經(jīng)過(guò)測(cè)試，這樣處理后速度基本沒(méi)有什么大的差異，這個(gè)方法也可以放棄。

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

? ? ?那么目前我想到的唯一的有可能對(duì)速度還有提升的就是定點(diǎn)化時(shí)不用32位的數(shù)據(jù)，適當(dāng)?shù)目紤]數(shù)據(jù)的范圍，如果能保證定點(diǎn)后的數(shù)據(jù)能在16位的有效范圍，那么還是有可能進(jìn)一步提高點(diǎn)速度的，畢竟這個(gè)時(shí)候可以使用SSE單指令一次性進(jìn)行8個(gè)整數(shù)的加減乘法了，這個(gè)有待于進(jìn)一步去測(cè)試。

　　六、特例優(yōu)化

　　在有些情況下甚至很多情況下，我們使用的Feature是其自身，這種情況下因?yàn)閿?shù)據(jù)的特殊性，我們可以做一些特殊處理，使得算法的速度更快。

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

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的灰度圖像計(jì)算用367ms左右，比上述都要快。

　　同樣的道理，我們也可以使用BCB技術(shù)來(lái)優(yōu)化，但是此時(shí)的BCB來(lái)的更簡(jiǎn)單。

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]++; // 計(jì)算每行第一個(gè)點(diǎn)的二維直方圖，直方圖的水平方向?yàn)镕eature坐標(biāo)，垂直方向?yàn)閂alue坐標(biāo) BCB[Value]--; // 此時(shí)的CutPoint初始化為-1，所以+方向的數(shù)據(jù)為0，所有的都在-方向 }}for (int X = 0; X < Width; X++){float BalanceWeight = 0;int IndexF = LinePS[X] * Level; // 中心點(diǎn)P的Value所對(duì)應(yīng)的那一行Feature權(quán)重起始索引for (int I = 0; I < Level; I++) // BCB[I]中保存的是以CutPoint為分界線，Feature為I時(shí)，分界線左側(cè)的所有Value[0-CutPoint]值的數(shù)量和分界線右側(cè)所有的Value（CutPoint， LevelV - 1]值數(shù)量的差異{BalanceWeight += BCB[I] * Weight[IndexF + I]; // 因?yàn)镕eature為固定值時(shí)，如果中心點(diǎn)固定，那么不管與Feature對(duì)應(yīng)的Value值時(shí)多少，Weight就是定值了。}if (BalanceWeight < 0) // 第一個(gè)點(diǎn)的BalanceWeight必然小于0{for (; BalanceWeight < 0 && CutPoint != Level - 1; CutPoint++){int Index = CutPoint + 1; // 新的直方圖的位置BCB[Index] += Histgram[Index] * 2; // 數(shù)量是同樣的道理BalanceWeight += 2 * Histgram[Index] * Weight[IndexF + Index];}}else if (BalanceWeight > 0) // 如果平衡值大于0，則向左移動(dòng)中間值{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++) // 即將移出的那一列數(shù)據(jù){int Value = Input[J * Stride + X - Radius];Histgram[Value]--;if (Value <= CutPoint) // 如果移出的那個(gè)值小于當(dāng)前的中值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) // 如果移出的那個(gè)值小于當(dāng)前的中值BCB[Value]++;elseBCB[Value]--;}}}} FreeMemory:if (Histgram != NULL) free(Histgram);if (BCB != NULL) free(BCB);return Status; }

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

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

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]++; // 計(jì)算每行第一個(gè)點(diǎn)的二維直方圖，直方圖的水平方向?yàn)镕eature坐標(biāo)，垂直方向?yàn)閂alue坐標(biāo) Balance--;}}for (int X = 0; X < Width; X++){ if (Balance < 0) // 第一個(gè)點(diǎn)的Balance必然小于0{for (; Balance < 0 && CutPoint != Level - 1; CutPoint++){ Balance += 2 * Histgram[CutPoint + 1];}}else if (Balance > 0) // 如果平衡值大于0，則向左移動(dòng)中間值{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++) // 即將移出的那一列數(shù)據(jù){int Value = Input[J * Stride + X - Radius];Histgram[Value]--;if (Value <= CutPoint) // 如果移出的那個(gè)值小于當(dāng)前的中值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) // 如果移出的那個(gè)值小于當(dāng)前的中值Balance++;elseBalance--;}}}} FreeMemory:if (Histgram != NULL) free(Histgram);return Status; }

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

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

? ? ?在論文的最后，講述了加權(quán)中值模糊的多個(gè)應(yīng)用場(chǎng)景，比如在光流、立體匹配、JPG瑕疵修復(fù)、藝術(shù)特效等等方面，我測(cè)試下幾個(gè)我能做的測(cè)試，確實(shí)有不錯(cuò)的效果，比如下面的JPG瑕疵修復(fù)。對(duì)簡(jiǎn)單的圖處理后確實(shí)蠻好的，如果在結(jié)合我之前研究的MLAA去鋸齒算法，恢復(fù)后的圖像質(zhì)量就更高了，如下所示：

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

? ?? ? ?? ? ?

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

　　　　　　　　　　　　原圖　　　　　　　　　　　　　　　　　　　　　　　　　　 ? ? ?　加權(quán)中值模糊（特征圖為原圖）　　　　　　　　　　　　　　　　　 MLAA后續(xù)處理后（邊緣更平滑）

? ? ? ? 另外，WMF的保邊特性感覺(jué)比其他的如導(dǎo)向?yàn)V波、雙邊濾波等等都要強(qiáng)烈的多，比如下圖： ?

? ?? ? ??

　　花朵的邊緣，下面的文字等等處理都還特別清晰，不像其他的保邊濾波器總有點(diǎn)模糊，這個(gè)特性也許用到一些增強(qiáng)上也會(huì)有很不錯(cuò)的效果。

　　按照上述文章的思路，我整理和編制一個(gè)簡(jiǎn)易的測(cè)試程序，用來(lái)論證論文和我博文中得一些數(shù)據(jù)，使用的VS2013編譯的，用C++做的DLL，C#做的UI測(cè)試界面，不依賴(lài)于任何其他第三方庫(kù)，目前只做了灰度圖的方案，因?yàn)椴噬脑?huà)也基本就是三個(gè)通道獨(dú)立寫(xiě)，可以通過(guò)拆分然后調(diào)用灰度的來(lái)實(shí)現(xiàn)。我也測(cè)試了下作者分享的VS工程，應(yīng)該比我提供的代碼速度稍微慢一點(diǎn)。

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

總結(jié)

以上是生活随笔為你收集整理的加权中值滤波的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。

中值

上一篇： 90后，一个即将成为程序员的我
下一篇： Oracle内存结构详解（三）--管理O