中值濾波(中值濾波計算)

對于數(shù)字圖像的去噪，前邊我們講了均值濾波算法與高斯濾波算法，此外很常見的還有中值濾波算法，這些濾波算法都屬于空間濾波，即對于每一個像素點(diǎn)，都選取其周圍矩形區(qū)域中的像素點(diǎn)來計算濾波值。最近在項目中要使用到中值濾波，發(fā)現(xiàn)如果調(diào)用Opencv的medianBlur函數(shù)來實(shí)現(xiàn)中值濾波，窗口為3*3或者5*5時耗時為幾毫秒，當(dāng)窗口達(dá)到7*7或者9*9以上，耗時將增加至幾十毫秒，這很影響實(shí)時性，所以自己基于C++、Opencv和CUDA實(shí)現(xiàn)了一個中值濾波函數(shù)，發(fā)現(xiàn)當(dāng)窗口達(dá)到7*7以上時，比Opencv快了不少，這還是很有成就感的。所以本文的主要內(nèi)容是講中值濾波的原理、實(shí)現(xiàn)和優(yōu)化。

1. 中值濾波的適用場景。中值濾波對鹽椒噪聲的去噪效果相當(dāng)好，其效果比均值濾波、高斯濾波要好得多。這里的鹽椒噪聲，我們也可以理解為一些隨機(jī)分布的、比較孤立的噪點(diǎn)，如下圖所示。在下面的內(nèi)容中，我們將使用不同的濾波算法分別對這張布滿鹽椒噪聲的圖像進(jìn)行去噪，并比較去噪結(jié)果。

噪聲圖

2. 中值濾波原理。其原理非常簡單，從字面就可以理解，所謂中值，就是對鄰域矩形窗口所有點(diǎn)的像素值進(jìn)行排序（從大到小或從小到大都o(jì)k），然后取排在最中間的像素值作為當(dāng)前待濾波點(diǎn)的濾波值。

假設(shè)取3*3窗口進(jìn)行中值濾波，如下圖所示，P4為當(dāng)前待濾波點(diǎn)，取其周圍矩形區(qū)域中9個點(diǎn)的像素值進(jìn)行排序，假設(shè)排序結(jié)果為：P0<P5<P7<P2<P3<P1<P8<P4<P6，那么則取P3的像素值作為P4的濾波值。

3. 查詢中值的算法選擇。查找中值的效率，是整個中值濾波算法效率的關(guān)鍵。所以為了盡可能地減少計算耗時，我們必須使用一種高效的查找中值算法。

(1) 冒泡排序。對矩形窗口內(nèi)所有點(diǎn)的像素值進(jìn)行從小到大或從大到小的冒泡排序，然后取排在最中間的像素值作為濾波值。這是最直接的中值查找算法，也是最耗時的算法。

冒泡排序的原理很簡單，就是一輪一輪地逐個比較，假設(shè)有9個像素值P0~P8，現(xiàn)使用冒泡排序?qū)ζ溥M(jìn)行從大到小地排序，排序過程分為8個步驟：

步驟一：分別依次比較P0與Pi（1 ≤ i ≤ 8）的大小。

如果任意Pi大于P0，都交換Pi與P0的值，假設(shè)P3>P0，那么則交換P3與P0的值：

步驟二：分別依次比較P1與Pi（2 ≤ i ≤ 8）的大小。如果任意Pi大于P1，都交換Pi與P1的值。

步驟三：分別依次比較P2與Pi（3 ≤ i ≤ 8）的大小。如果任意Pi大于P2，都交換Pi與P2的值。

.

.

.

步驟八：比較P7與P8的大小。如果P8大于P7，則交換P8與P7的值。

(2) 快速排序。

快速排序的核心思想：首先把所有數(shù)據(jù)看成一類，選擇一個界線值P0（通常把待排序數(shù)據(jù)的第一個值P0作為界線值），然后以P0為分界線，把其余大于P0的數(shù)據(jù)放到P0的左側(cè)，小于等于P0的數(shù)據(jù)放到P0的右側(cè)。此時P0左、右側(cè)的數(shù)據(jù)分別組成一個新的類，然后繼續(xù)對新類進(jìn)行以上同樣的操作，一直到所有類的數(shù)據(jù)個數(shù)都為1為止。

所以快速排序是一個裂變的過程：

(3) 分行求中值。針對于中值濾波算法，我們的目的是求取矩形窗口中所有像素值的中值，實(shí)際上只要求得中值就ok了，并不需要像上面一樣對窗口中所有像素值進(jìn)行排序。假設(shè)取3*3窗口，窗口內(nèi)的所有像素值為P0~P8，P0~P2為窗口的一行，P3~P5為窗口的一行，P6~P8為窗口的一行，那么我們可以按照以下步驟求取3*3窗口的中值：

步驟一：對P0~P2進(jìn)行排序，得到這三個數(shù)中的最小值、中值、最大值分別記為min0、mid0、max0。

步驟二：對P3~P5進(jìn)行排序，得到這三個數(shù)中的最小值、中值、最大值分別記為min1、mid1、max1。

步驟三：對P6~P8進(jìn)行排序，得到這三個數(shù)中的最小值、中值、最大值分別記為min2、mid2、max2。

步驟四：求min0、min1、min2的最大值max(min)。

步驟五：求mid0、mid1、mid2的中值mid(mid)。

步驟六：求max0、max1、max2的最小值min(max)。

步驟七：求max(min)、mid(mid)、min(max)的中值，即為P0~P8的中值。

這種算法雖然減少了不少計算量，但是當(dāng)窗口尺寸增加之后，多次求局部最小值、中值、最大值也是比較費(fèi)力的，整體耗時也就上來了。

(4) 使用窗口的灰度直方圖來求取中值。

主要分為兩個步驟：

步驟一：統(tǒng)計窗口內(nèi)的灰度直方圖。

步驟二：使用上一步得到的灰度直方圖，從第0級灰度開始計算累加直方圖。當(dāng)某一級灰度i的累加直方圖值達(dá)到窗口總像素的一半時，則判定i就是窗口內(nèi)所有像素值的中值。

圖像的灰度級是確定的，比如8位圖像，其灰度級（像素值）為0~255，那么首先統(tǒng)計窗口內(nèi)0~255像素值各自的點(diǎn)數(shù)，即為統(tǒng)計直方圖，然后從0開始往后計算各個灰度級的累加直方圖值，累加直方圖值即對于每一灰度級，窗口內(nèi)像素值小于等于該灰度級的點(diǎn)數(shù)。比如當(dāng)計算到灰度級50時，判定發(fā)現(xiàn)灰度級50的的累加直方圖達(dá)到窗口總像素的一半，則停止計算，取50作為中值。

綜合考慮，使用直方圖獲取中值，計算復(fù)雜度相對更低，因此我們使用該方法來實(shí)現(xiàn)中值濾波。

4. 代碼實(shí)現(xiàn)。

首先是CUDA的核函數(shù)代碼：

#define WIN_SIZE 3#define WIN_SIZE_2 (WIN_SIZE>>1)#defineD_SIZE(WIN_SIZE*WIN_SIZE)#define LEVEL 256 //0~255//定義紋理內(nèi)存，用于快速訪問源圖像texture<uchar, cudaTextureType2D, cudaReadModeElementType> tex_src;//源圖像使用紋理內(nèi)存綁定訪問__global__ void mediablur_ker(uchar *out, int row, int col){ int x = blockIdx.x * blockDim.x + threadIdx.x; //col int y = blockIdx.y * blockDim.y + threadIdx.y; //row if (x < col && y < row){ uchar hist[LEVEL] = {0}; for (int i = 0; i < WIN_SIZE; i++) //統(tǒng)計窗口內(nèi)的灰度直方圖 { for (int j = 0; j < WIN_SIZE; j++) { uchar idx = tex2D(tex_src, (x+j), (y+i)); hist[idx]++; } } uchar mid = 0; uchar n = D_SIZE / 2 + 1; for (int i = 0; i < LEVEL; i++) //計算累加直方圖 { mid += hist[i]; if (mid >= n) //當(dāng)某一級灰度i的累加直方圖值達(dá)到窗口總像素的一半時，則判定i就是窗口內(nèi)所有像素值的中值 { out[y*col + x] = (uchar)i; break; } } }}

其次是實(shí)現(xiàn)函數(shù)代碼：

void mediablur_fiter_cuda(Mat src, Mat &dst){ Mat src_board; //擴(kuò)充邊緣 copyMakeBorder(src, src_board, WIN_SIZE_2, WIN_SIZE_2, WIN_SIZE_2, WIN_SIZE_2, BORDER_REFLECT); //擴(kuò)充邊緣 const int row0 = src.rows; const int col0 = src.cols; const int img_size0 = row0*col0; const int row = src_board.rows; const int col = src_board.cols; const int img_size = row*col; ////////////////////////////////////////////////////////////////////////////////////////////////////////////// uchar *dst_cuda; cudaMalloc((void**)&dst_cuda, img_size0); ////////////////////////////////////////////////////////////////////////////////////////////////////////////// cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<uchar>();//聲明數(shù)據(jù)類型 cudaArray *cuArray_src; cudaMallocArray(&cuArray_src, &channelDesc, col, row); //分配大小為col*row的CUDA數(shù)組 tex_src.addressMode[0] = cudaAddressModeWrap;//尋址方式 tex_src.addressMode[1] = cudaAddressModeWrap;//尋址方式 如果是三維數(shù)組則設(shè)置texRef.addressMode[2] tex_src.normalized = fal;//是否對紋理坐標(biāo)歸一化 tex_src.filterMode = cudaFilterModePoint;//紋理的濾波模式：最近點(diǎn)取樣和線性濾波 cudaFilterModeLinear cudaBindTextureToArray(&tex_src, cuArray_src, &channelDesc); //紋理綁定，CUDA數(shù)組和紋理參考的連接 cudaMemcpyToArray(cuArray_src, 0, 0, src_board.data, img_size, cudaMemcpyHostToDevice); ////////////////////////////////////////////////////////////////////////////////////////////////////////////// dim3 Block_G(16, 16); dim3 Grid_G((col0 + 15) / 16, (row0 + 15) / 16); //調(diào)用核函數(shù) mediablur_ker << <Grid_G, Block_G >> >(dst_cuda, row0, col0); ////////////////////////////////////////////////////////////////////////////////////////////////////////////// dst = Mat::zeros(src.size(), CV_8UC1); cudaMemcpy(dst.data, dst_cuda, img_size0, cudaMemcpyDeviceToHost); //cout << dst << endl; ////////////////////////////////////////////////////////////////////////////////////////////////////////////// cudaFree(dst_cuda); cudaFreeArray(cuArray_src); cudaUnbindTexture(tex_src); }

5. 濾波結(jié)果對比。

(1) 統(tǒng)一取3*3窗口，分別調(diào)用以上實(shí)現(xiàn)的中值濾波函數(shù)，和Opencv的高斯濾波函數(shù)、均值濾波函數(shù)對本文開頭的噪聲圖進(jìn)行濾波去噪，得到的結(jié)果如下圖所示。可以看到，中值濾波對鹽椒噪聲的去噪效果，相比于高斯濾波和均值濾波來說，那是相當(dāng)出色。

均值濾波結(jié)果

高斯濾波結(jié)果

中值濾波結(jié)果

(2) 比較以上實(shí)現(xiàn)的中值濾波與Opencv實(shí)現(xiàn)的中值濾波函數(shù)的耗時。如下表所示，可以看到，當(dāng)窗口在7*7以上時，本文的優(yōu)化加速效果就很明顯了。

3*3

5*5

7*7

9*9

本文實(shí)現(xiàn)/ms

2.69

2.985

3.466

3.930

Opencv函數(shù)/ms

0.475

2.243

28.093

32.197