中值滤波能够有效去除图像中的异常点,具有去除图像噪声的作用。传统中值滤波的算法一般都是在图像中建立窗口,然后对窗口内的所有像素值进行排序,选择排序后的中间值作为窗口中心像素滤波后的值。由于这个做法在每个像素点处都要建立窗口并排序,非常耗时,尤其是有大量的冗余计算。如下图:
黄色区域+中间粉色区域是第一个像素为中心建立的滤波窗口,粉色区域+右边蓝色区域为同一行第二个像素为中心建立的滤波窗口。传统做法对每一个窗口内所有像素排序,那么中间粉色区域的像素就会被重复的做排序运算,存在大量冗余计算。如果窗口移动一个像素的时候只用考虑去除左边一列内(黄色区域)的像素,并且加上右边一列(蓝色区域)的像素,运算会大大简化。这样的操作可以使用直方图来实现。
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中值滤波结果:
直方图快速滤波结果:
完整代码如下:(博主偷懒没有提前做边界扩展,而是直接保留了四个边界的像素值,边界扩展也很容易实现,不再赘述)
优化效果:在我年久失修的I3笔记本上跑512*512大小的图像,OpenCV中的Mediablur()需要约200ms,优化后的fastMediaBlur()只需要约70ms。不同图像,不同平台时间会有差异,但是速度确实提升约2-3倍。
- #include <opencv2\opencv.hpp>
- #include <iostream>
- #include <string>
- using namespace cv;
- using namespace std;
- //计算亮度中值和灰度<=中值的像素点个数
- void CalculateImage_MedianGray_PixelCount(const Mat &histogram, int pixelCount, int &medianValue, int &pixleCountLowerMedian)
- {
- float *data = (float *)histogram.data;//直方图
- int sum = 0;
- for (int i = 0; i <= 255; ++i)
- {
- //
- sum += data[i];
- if (2 * sum>pixelCount)
- {
- medianValue = i;
- pixleCountLowerMedian = sum;
- break;
- }
- }
- }
- void fastMedianBlur(const Mat &srcImg, Mat &dstImg, int diameter)
- {
- int radius = (diameter - 1) / 2;
- int imgW = srcImg.cols;
- int imgH = srcImg.rows;
- int channels = srcImg.channels();
- dstImg = Mat::zeros(imgH, imgW, CV_8UC1);
- int windowSize = diameter*diameter;
- Mat windowImg = Mat::zeros(diameter, diameter, CV_8UC1);
- //直方图
- Mat histogram;
- int histogramSize = 256;//灰度等级
- int thresholdValue = windowSize / 2 + 1;//step1.设置阈值(步骤参考:图像的高效编程要点之四)
- //待处理图像像素
- uchar *pDstData = dstImg.data + imgW*radius + radius;
- //整个图像中窗口位置指针
- uchar *pSrcData = srcImg.data;
- //逐行遍历
- for (int i = radius; i <= imgH - 1 - radius; i++)
- {
- //列指针
- uchar *pColDstData = pDstData;
- uchar *pColSrcData = pSrcData;
- //单个窗口指针
- uchar *pWindowData = windowImg.data;
- uchar *pRowSrcData = pColSrcData;
- //从源图中提取窗口图像
- for (int winy = 0; winy <= diameter - 1; winy++)
- {
- for (int winx = 0; winx <= diameter - 1; winx++)
- {
- pWindowData[winx] = pRowSrcData[winx];
- }
- pRowSrcData += imgW;
- pWindowData += diameter;
- }
- //求直方图,确定中值,并记录亮度<=中值的像素点个数
- calcHist(&windowImg,
- 1,//Mat的个数
- 0,//用来计算直方图的通道索引,通道索引依次排开
- Mat(),//Mat()返回一个空值,表示不用mask,
- histogram, //直方图
- 1, //直方图的维数,如果计算2个直方图,就为2
- &histogramSize, //直方图的等级数(如灰度等级),也就是每列的行数
- 0//分量的变化范围
- );
- //求亮度中值和<=中值的像素点个数
- int medianValue, pixelCountLowerMedian;
- CalculateImage_MedianGray_PixelCount(histogram, windowSize, medianValue, pixelCountLowerMedian);
- ////////////滑动窗口操作结束///////////////////////
- //滤波
- pColDstData[0] = (uchar)medianValue;
- //处理同一行下一个像素
- pColDstData++;
- pColSrcData++;
- for (int j = radius + 1; j <= imgW - radius - 1; j++)
- {
- //维护滑动窗口直方图
- //删除左侧
- uchar *pWinLeftData = pColSrcData - 1;
- float *pHistData = (float*)histogram.data;
- for (int winy = 0; winy < diameter; winy++)
- {
- uchar grayValue = pWinLeftData[0];
- pHistData[grayValue] -= 1.0;
- if (grayValue <= medianValue)
- {
- pixelCountLowerMedian--;
- }
- pWinLeftData += imgW;
- }
- //增加右侧
- uchar *pWinRightData = pColSrcData + diameter - 1;
- for (int winy = 0; winy < diameter; winy++)
- {
- uchar grayValue = pWinRightData[0];
- pHistData[grayValue] += 1.0;
- if (grayValue <= medianValue)
- {
- pixelCountLowerMedian++;
- }
- pWinRightData += imgW;
- }
- //计算新的中值
- if (pixelCountLowerMedian > thresholdValue)
- {
- while (1)
- {
- pixelCountLowerMedian -= pHistData[medianValue];
- medianValue--;
- if (pixelCountLowerMedian <= thresholdValue)
- {
- break;
- }
- }
- }
- else
- {
- while (pixelCountLowerMedian < thresholdValue)
- {
- medianValue++;
- pixelCountLowerMedian += pHistData[medianValue];
- }
- }
- pColDstData[0] = medianValue;
- //下一个像素
- pColDstData++;
- pColSrcData++;
- }
- //移动至下一行
- pDstData += imgW;
- pSrcData += imgW;
- }
- //边界直接赋原始值,不做滤波处理
- pSrcData = srcImg.data;
- pDstData = dstImg.data;
- //上下边界
- for (int j = 0; j < imgW; j++)
- {
- for (int i = 0; i < radius; i++)
- {
- int idxTop = i*imgW + j;
- pDstData[idxTop] = pSrcData[idxTop];
- int idxBot = (imgH - i - 1)*imgW + j;
- pDstData[idxBot] = pSrcData[idxBot];
- }
- }
- //左右边界
- for (int i = radius; i < imgH - radius - 1; i++)
- {
- for (int j = 0; j < radius; j++)
- {
- int idxLeft = i*imgW + j;
- pDstData[idxLeft] = pSrcData[idxLeft];
- int idxRight = i*imgW + imgW - j-1;
- pDstData[idxRight] = pSrcData[idxRight];
- }
- }
- }
- void main()
- {
- string imgPath = "data/";
- Mat srcImg = imread(imgPath + "moon.jpg", 0);
- Mat dstImg;
- double t0 = cv::getTickCount();
- fastMedianBlur(srcImg, dstImg, 5);
- //cv::medianBlur(srcImg, dstImg, 5); //OpenCV
- double t1 = cv::getTickCount();
- cout << "time=" << (t1 - t0) / cv::getTickFrequency() << endl;
- imwrite("data/test/srcImg.bmp", srcImg);
- imwrite("data/test/myFilter.bmp", dstImg);
- }
来源: http://www.cnblogs.com/riddick/p/7989871.html