中值滤波能够有效去除图像中的异常点,具有去除图像噪声的作用。传统中值滤波的算法一般都是在图像中建立窗口,然后对窗口内的所有像素值进行排序,选择排序后的中间值作为窗口中心像素滤波后的值。由于这个做法在每个像素点处都要建立窗口并排序,非常耗时,尤其是有大量的冗余计算。如下图:
黄色区域+中间粉色区域是第一个像素为中心建立的滤波窗口,粉色区域+右边蓝色区域为同一行第二个像素为中心建立的滤波窗口。传统做法对每一个窗口内所有像素排序,那么中间粉色区域的像素就会被重复的做排序运算,存在大量冗余计算。如果窗口移动一个像素的时候只用考虑去除左边一列内(黄色区域)的像素,并且加上右边一列(蓝色区域)的像素,运算会大大简化。这样的操作可以使用直方图来实现。
1.读取图像I,并且设定滤波窗口大小(winX*winY),一般winX=winY,奇数。
2.设定中值滤波直方图中的阈值,Thresh=(winX*winY)/2 +1;
3.如果要考虑边界情况,可以先对原图像进行扩展,左、右边界分别扩展winX/2个像素,上下边界分别扩展winY/2个像素。
4.逐行遍历图像像素,以第一行为例:先取第一行第一个要处理的像素(窗口中心像素),建立滤波窗口,提取窗口内所有像素值(N=winX*winY个),获取N个像素的直方图Hist。从左到右累加直方图中每个灰度层级下像素点个数,记为sumCnt,直到sumCnt>=Thresh,这时的灰度值就是当前窗口内所有像素值的中值MediaValue。将MediaValue值赋值给窗口中心像素,表明第一个像素中值滤波完成。
5.此时窗口需要向右移动一个像素,开始滤波第二个像素,并且更新直方图。以第二个像素为窗口中心建立滤波窗口,从前一个窗口的灰度直方图Hist中减去窗口中最左侧的一列像素值的灰度个数,然后加上窗口最右侧一列像素值的灰度个数。完成直方图的更新。
6.直方图更新后,sumCnt值有三种变化可能:(1)减小(2)维持不变(3)增大。这三种情况与减去与加入的像素值灰度有关。此时为了求得新的中值,需要不断调整sumCnt与Thresh之间的关系。
(1)如果sumCnt值小于Thresh:说明中值在直方图当前灰度层级的右边,sumCnt就依次向右加上一个灰度层级中灰度值个数,直到满足sumCnt>=Thresh为止。记录此时的灰度层级代表的灰度值,更新MediaValue,作为第二个像素的滤波后的值。
(2)维持不变:说明MediaValue值不变,直接作为第二个像素滤波后的值。
(3)如果sumCnt值大于Thresh:说明中值在直方图当前灰度层级的左边,sumCnt就依次向左减去一个灰度层级中灰度值个数,直到满足sumCnt<=Thresh为止。记录此时的灰度层级代表的灰度值,更新MediaValue值,作为第二个像素的滤波后的值。
7.窗口逐行依次滑动,求得整幅图像的中值滤波结果。
二、 滤波结果
以下图手机拍摄的moon.jpg为例:
OpenCV中值滤波结果:
直方图快速滤波结果:
完整代码(两种实现,原理一样)如下:(博主偷懒没有提前做边界扩展,而是直接保留了四个边界的像素值,边界扩展也很容易实现,不再赘述)
Code01:
1 #include <opencv2\opencv.hpp>
2 #include <iostream>
3 #include <string>
4
5 using namespace cv;
6 using namespace std;
7
8 //计算亮度中值和灰度<=中值的像素点个数
9 int calMediaValue(const int histogram[], int thresh)
10 {
11 int sum = 0;
12 for (int i = 0; i < 256; ++i)
13 {
14 sum += histogram[i];
15 if (sum>= thresh)
16 {
17 return i;
18 }
19 }
20 return 255;
21 }
22
23 void fastMedianBlur(const Mat &srcImg, Mat &dstImg, int diameter)
24 {
25 int radius = (diameter - 1) / 2;
26 int imgW = srcImg.cols;
27 int imgH = srcImg.rows;
28 int channels = srcImg.channels();
29 dstImg = Mat::zeros(imgH, imgW, CV_8UC1);
30 int windowSize = diameter*diameter;
31
32 //直方图
33 int Hist[256]={0};
34 int histogramSize = 256;//灰度等级
35 int thresholdValue = windowSize / 2 + 1;
36
37 uchar *pSrcData=srcImg.data;
38 uchar *pDstData=dstImg.data;
39
40 int right=imgW-radius;
41 int bot=imgH-radius;
42 for (int i=radius; i<right; i++)
43 {
44 for (int j=radius; j<bot; j++)
45 {
46 //每一行第一个待滤波像素建立直方图
47 if(j==radius)
48 {
49 memset(Hist, 0, histogramSize*sizeof(int));
50 for (int y=i-radius; y<=i+radius; y++)
51 {
52 for (int x=j-radius; x<=j+radius; x++)
53 {
54 uchar value=pSrcData[ y*imgW+x];
55 Hist[value]++;
56 }
57 }
58 }
59 else//更新直方图
60 {
61 int left=j-radius-1;
62 int right=j+radius;
63 for (int y=i-radius; y<=i+radius; y++)
64 {
65 //减去左边一列
66 int leftIdx=y*imgW+left;
67 uchar leftValue=pSrcData[leftIdx];
68 Hist[leftValue]--;
69
70 //加上右边一列
71 int rightIdx=y*imgW+right;
72 uchar rightValue=pSrcData[rightIdx];
73 Hist[rightValue]++;
74 }
75 }
76
77 //直方图求中值
78 uchar filterValue=calMediaValue(Hist, thresholdValue);
79 pDstData[i*imgW+j]=filterValue;
80 }
81 }
82
83 //边界直接赋原始值,不做滤波处理
84 pSrcData = srcImg.data;
85 pDstData = dstImg.data;
86 //上下边界
87 for (int j = 0; j < imgW; j++)
88 {
89 for (int i = 0; i < radius; i++)
90 {
91 int idxTop = i*imgW + j;
92 pDstData[idxTop] = pSrcData[idxTop];
93 int idxBot = (imgH - i - 1)*imgW + j;
94 pDstData[idxBot] = pSrcData[idxBot];
95 }
96 }
97 //左右边界
98 for (int i = radius; i < imgH - radius - 1; i++)
99 {
100 for (int j = 0; j < radius; j++)
101 {
102 int idxLeft = i*imgW + j;
103 pDstData[idxLeft] = pSrcData[idxLeft];
104 int idxRight = i*imgW + imgW - j-1;
105 pDstData[idxRight] = pSrcData[idxRight];
106 }
107 }
108 }
109
110
111 void main()
112 {
113 string imgPath = "data/";
114 Mat srcImg = imread(imgPath + "moon.jpg", 0);
115 Mat dstImg;
116 double t0 = cv::getTickCount();
117 fastMedianBlur(srcImg, dstImg, 5);
118 double t1 = cv::getTickCount();
119 cout << "time=" << (t1 - t0) / cv::getTickFrequency() << endl;
120
121 imwrite("data/test/srcImg.jpg", srcImg);
122 imwrite("data/test/myFilter.jpg", dstImg);
123 }
Code02:
1 #include <opencv2\opencv.hpp>
2 #include <iostream>
3 #include <string>
4
5 using namespace cv;
6 using namespace std;
7
8 //计算亮度中值和灰度<=中值的像素点个数
9 void CalculateImage_MedianGray_PixelCount(const Mat &histogram, int pixelCount, int &medianValue, int &pixleCountLowerMedian)
10 {
11 float *data = (float *)histogram.data;//直方图
12 int sum = 0;
13 for (int i = 0; i <= 255; ++i)
14 {
15 //
16 sum += data[i];
17 if (2 * sum>pixelCount)
18 {
19 medianValue = i;
20 pixleCountLowerMedian = sum;
21 break;
22 }
23 }
24 }
25
26 void fastMedianBlur(const Mat &srcImg, Mat &dstImg, int diameter)
27 {
28 int radius = (diameter - 1) / 2;
29 int imgW = srcImg.cols;
30 int imgH = srcImg.rows;
31 int channels = srcImg.channels();
32 dstImg = Mat::zeros(imgH, imgW, CV_8UC1);
33 int windowSize = diameter*diameter;
34 Mat windowImg = Mat::zeros(diameter, diameter, CV_8UC1);
35
36 //直方图
37 Mat histogram;
38 int histogramSize = 256;//灰度等级
39 int thresholdValue = windowSize / 2 + 1;//step1.设置阈值(步骤参考:图像的高效编程要点之四)
40
41 //待处理图像像素
42 uchar *pDstData = dstImg.data + imgW*radius + radius;
43 //整个图像中窗口位置指针
44 uchar *pSrcData = srcImg.data;
45
46 //逐行遍历
47 for (int i = radius; i <= imgH - 1 - radius; i++)
48 {
49 //列指针
50 uchar *pColDstData = pDstData;
51 uchar *pColSrcData = pSrcData;
52
53 //单个窗口指针
54 uchar *pWindowData = windowImg.data;
55 uchar *pRowSrcData = pColSrcData;
56 //从源图中提取窗口图像
57 for (int winy = 0; winy <= diameter - 1; winy++)
58 {
59 for (int winx = 0; winx <= diameter - 1; winx++)
60 {
61 pWindowData[winx] = pRowSrcData[winx];
62 }
63 pRowSrcData += imgW;
64 pWindowData += diameter;
65 }
66
67 //求直方图,确定中值,并记录亮度<=中值的像素点个数
68 calcHist(&windowImg,
69 1,//Mat的个数
70 0,//用来计算直方图的通道索引,通道索引依次排开
71 Mat(),//Mat()返回一个空值,表示不用mask,
72 histogram, //直方图
73 1, //直方图的维数,如果计算2个直方图,就为2
74 &histogramSize, //直方图的等级数(如灰度等级),也就是每列的行数
75 0//分量的变化范围
76 );
77
78 //求亮度中值和<=中值的像素点个数
79 int medianValue, pixelCountLowerMedian;
80 CalculateImage_MedianGray_PixelCount(histogram, windowSize, medianValue, pixelCountLowerMedian);
81 ////////////滑动窗口操作结束///////////////////////
82
83 //滤波
84 pColDstData[0] = (uchar)medianValue;
85
86 //处理同一行下一个像素
87 pColDstData++;
88 pColSrcData++;
89 for (int j = radius + 1; j <= imgW - radius - 1; j++)
90 {
91 //维护滑动窗口直方图
92 //删除左侧
93 uchar *pWinLeftData = pColSrcData - 1;
94 float *pHistData = (float*)histogram.data;
95 for (int winy = 0; winy < diameter; winy++)
96 {
97 uchar grayValue = pWinLeftData[0];
98 pHistData[grayValue] -= 1.0;
99 if (grayValue <= medianValue)
100 {
101 pixelCountLowerMedian--;
102 }
103 pWinLeftData += imgW;
104 }
105
106 //增加右侧
107 uchar *pWinRightData = pColSrcData + diameter - 1;
108 for (int winy = 0; winy < diameter; winy++)
109 {
110 uchar grayValue = pWinRightData[0];
111 pHistData[grayValue] += 1.0;
112 if (grayValue <= medianValue)
113 {
114 pixelCountLowerMedian++;
115 }
116 pWinRightData += imgW;
117 }
118 //计算新的中值
119 if (pixelCountLowerMedian > thresholdValue)
120 {
121 while (1)
122 {
123 pixelCountLowerMedian -= pHistData[medianValue];
124 medianValue--;
125 if (pixelCountLowerMedian <= thresholdValue)
126 {
127 break;
128 }
129 }
130 }
131 else
132 {
133 while (pixelCountLowerMedian < thresholdValue)
134 {
135 medianValue++;
136 pixelCountLowerMedian += pHistData[medianValue];
137
138 }
139 }
140
141 pColDstData[0] = medianValue;
142 //下一个像素
143 pColDstData++;
144 pColSrcData++;
145 }
146 //移动至下一行
147 pDstData += imgW;
148 pSrcData += imgW;
149 }
150
151 //边界直接赋原始值,不做滤波处理
152 pSrcData = srcImg.data;
153 pDstData = dstImg.data;
154 //上下边界
155 for (int j = 0; j < imgW; j++)
156 {
157 for (int i = 0; i < radius; i++)
158 {
159 int idxTop = i*imgW + j;
160 pDstData[idxTop] = pSrcData[idxTop];
161 int idxBot = (imgH - i - 1)*imgW + j;
162 pDstData[idxBot] = pSrcData[idxBot];
163 }
164 }
165 //左右边界
166 for (int i = radius; i < imgH - radius - 1; i++)
167 {
168 for (int j = 0; j < radius; j++)
169 {
170 int idxLeft = i*imgW + j;
171 pDstData[idxLeft] = pSrcData[idxLeft];
172 int idxRight = i*imgW + imgW - j-1;
173 pDstData[idxRight] = pSrcData[idxRight];
174 }
175 }
176 }
177
178
179 void main()
180 {
181 string imgPath = "data/";
182 Mat srcImg = imread(imgPath + "moon.jpg", 0);
183 Mat dstImg;
184 double t0 = cv::getTickCount();
185 fastMedianBlur(srcImg, dstImg, 5);
186 //cv::medianBlur(srcImg, dstImg, 5); //OpenCV
187 double t1 = cv::getTickCount();
188 cout << "time=" << (t1 - t0) / cv::getTickFrequency() << endl;
189
190 imwrite("data/test/srcImg.bmp", srcImg);
191 imwrite("data/test/myFilter.bmp", dstImg);
192 }