【影象處理】高效的中值濾波(bug 已修復）
之前的bug：
當灰度為255或者0時，出現灰度溢位的bug，導致灰度黑白顛倒，現已修復，並重新將函式改為無返回值型別，原有的帶有返回影象的函式不規範，容易忘記釋放空間。新的函式再最後面
經測試，我的程式計算速度比OpenCV耗時多多了，，，~~(>_<)~~ ，，真心比不起OpenCV，因為這貨用intel 的IPP庫了，誰讓人家是親兒子呢，，，，
（1）中值濾波
由於中值濾波是一個非線性濾波器，中值的求取需要排序，傳統的中值濾波演算法中無非是在排序上做文章，爭取做到最優的快速排序，但是本文所實現的方法不同，該方法充分利用了直方圖。
從程式設計的觀點看，直方圖是一種很有效的資料結構，所佔記憶體空間很少，又能反映出影象中的灰度分佈和目標特性等等，且直方圖本身就是有序的。基於直方圖可以很容易、很高效地得到影象中亮度、對比度、最大亮度、最小亮度及亮度中值。
本文章旨在實現高效的中值濾波（對於這個大眾化的影象增強方法沒有深究，也不知該方法到底源自哪一片文章，鄙人知識有幸找到了一份老師給的PDF文件，上面介紹了該方法，就抽時間實現了一下，效果和OpenCV裡實現的效果是一樣的，至於速度，只會快，不會慢）
（2）原理
如上圖所示，在進行橫向滑動視窗濾波時，視窗中的畫素僅僅是丟掉了左側一列，增加了右側一列資料，如果丟掉中間重疊的這一部分資料，到下個視窗再重新定址和讀取資料，無疑是計算的沉重負擔，所以，該演算法的核心思想就是充分利用重疊部分，使用直方圖來計算中值，不需要排序演算法，快，且高效。
（3）演算法流程
假設視窗尺寸為 N × M
- 設定門限
th=N∗M2+1 - 將視窗移動到一個新行的開始，建立視窗畫素的直方圖，通過直方圖確定中值 med，記下亮度小於或等於med 的畫素數目到 mNum 中
對於最左列的每個畫素，設其灰度為
gl ，執行：（即去掉左側）
H[gl]= H[gl]−1若gl≤med，則mNum=mNum−1 將視窗右移一列，對於最右列的每個畫素，設其灰度為
gr ，執行：（即添上右側）
H[gr]=H[gr]+1 若gr≤med，則mNum=mNum+1 若
mNum>th ,則轉 6；否則(mNum≤th )：
while(mNum<th) {med=med+1;mNum=mNum+H[med];} 轉 7 重複
mNum=mNum−H[med] med=med−1 直到
mNum≤th - 如果視窗的右側不是影象的右邊界，則轉 3
- 若果視窗的底行不是圖象的下邊界，則轉 2
（4）實現
/**
* @brief 其中boundary1(src, width, height, mask_width, mask_height, &width1, &height1);是我自己實現的根據滑動視窗尺寸對原影象的邊界進行擴充套件的函式，擴充套件方式採用複製最臨近的畫素的方法進行，此處就不貼出來了。
getHist_local(pCurOrg, hist, width1, mask_width, mask_height)是實現的區域性直方圖函式，pCurOrg為要統計直方圖的區域的左上角畫素地址，width1是原影象的行步長。
/**
* @brief Get median grayscale
* @param histogram
* @param imgSize
* @return medianGray
*/
unsigned char getMedianGry(int *histogram, int imgSize)
{
int sum;
unsigned char medGry;
for(sum = 0, medGry = 0; medGry < 256; medGry++)
{
sum += histogram[medGry];
if(sum * 2 > imgSize)break;
}
return medGry;
}該函式實現通過直方圖尋找中值的功能，在上一個程式中有呼叫。
（5）結果
原影象（女神圖）
中值濾波後的影象（只抽取並處理了灰度圖）
處理結果與photoshop與OpenCV以及Matlab中基本上都是一致的。
改進和修復bug的版本
/**
* @brief Fast median filter
* @param src
* @param width
* @param height
* @param mask_width
* @param mask_height
*/
void medianBlur(unsigned char* src, int width, int height, int mask_width, int mask_height)
{
int hist[256];
int i, j, k, nMin;
int width1 = width + mask_width - 1;
int height1 = height + mask_height - 1;
int th = mask_width * mask_height / 2 + 1;
unsigned char *pCurOrg, *pCurObj, *pl, *pr;
unsigned char med;
unsigned char *src1 = Malloc(unsigned char, width1 * height1);
//獲得帶邊界擴充套件的影象副本
boundaryExtention(src, src1, width, height, mask_width, mask_height, width1, height1);
//清空原影象資料
memset(src, 0, width * height * sizeof(unsigned char));
for(i = 0; i < height; i++)
{
nMin = 0;
pCurObj = src + i * width;
pCurOrg = src1 + i * width1;
getHistLocal(pCurOrg, hist, width1, mask_width, mask_height);
med = getMedianGry(hist, mask_width * mask_height);
for(j = 0; j <= med; j++)
nMin += hist[j];
(*pCurObj++) = med;
for(k = 0; k < width - 1; k++, pCurOrg++, pCurObj++)
{
for(j = 0, pl = pCurOrg, pr = pCurOrg + mask_width; j < mask_height; j++, pl += width1, pr += width1)
{
if(*pl <= med ) nMin -= 1;
if(*pr <= med ) nMin += 1;
hist[*pl] -= 1;
hist[*pr] += 1;
}
while(nMin > th)
{
if(med==0)
break;
else
{
nMin -= hist[med];
med--;
}
}
while(nMin < th)
{
if(med==255)
break;
else
{
med++;
nMin += hist[med];
}
}
*pCurObj = med;
}
}
free(src1);
return;
}