中值滤波能够有效去除图像中的异常点，具有去除图像噪声的作用。传统中值滤波的算法一般都是在图像中建立窗口，然后对窗口内的所有像素值进行排序，选择排序后的中间值作为窗口中心像素滤波后的值。由于这个做法在每个像素点处都要建立窗口并排序，非常耗时，尤其是有大量的冗余计算。如下图：
黄色区域+中间粉色区域是第一个像素为中心建立的滤波窗口，粉色区域+右边蓝色区域为同一行第二个像素为中心建立的滤波窗口。传统做法对每一个窗口内所有像素排序，那么中间粉色区域的像素就会被重复的做排序运算，存在大量冗余计算。如果窗口移动一个像素的时候只用考虑去除左边一列内（黄色区域）的像素，并且加上右边一列（蓝色区域）的像素，运算会大大简化。这样的操作可以使用直方图来实现。
一、直方图实现快速中值滤波算法流程：
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.窗口逐行依次滑动，求得整幅图像的中值滤波结果。
具体代码如下：
#include <opencv2\opencv.hpp>
#include <iostream>
#include <string>
using namespace cv;
using namespace std;
//计算亮度中值和灰度<=中值的像素点个数
int calMediaValue(const int histogram[], int thresh)
{
int sum = 0;
for (int i = 0; i < 256; ++i)
{
sum += histogram[i];
if (sum>= thresh)
{
return i;
}
}
return 255;
}
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;
//直方图
int Hist[256]={0};
int histogramSize = 256;//灰度等级
int thresholdValue = windowSize / 2 + 1;
uchar *pSrcData=srcImg.data;
uchar *pDstData=dstImg.data;
int right=imgW-radius;
int bot=imgH-radius;
for (int i=radius; i<right; i++)
{
for (int j=radius; j<bot; j++)
{
//每一行第一个待滤波像素建立直方图
if(j==radius)
{
memset(Hist, 0, histogramSize*sizeof(int));
for (int y=i-radius; y<=i+radius; y++)
{
for (int x=j-radius; x<=j+radius; x++)
{
uchar value=pSrcData[ y*imgW+x];
Hist[value]++;
}
}
}
else//更新直方图
{
int left=j-radius-1;
int right=j+radius;
for (int y=i-radius; y<=i+radius; y++)
{
//减去左边一列
int leftIdx=y*imgW+left;
uchar leftValue=pSrcData[leftIdx];
Hist[leftValue]--;
//加上右边一列
int rightIdx=y*imgW+right;
uchar rightValue=pSrcData[rightIdx];
Hist[rightValue]++;
}
}
//直方图求中值
uchar filterValue=calMediaValue(Hist, thresholdValue);
pDstData[i*imgW+j]=filterValue;
}
}
//边界直接赋原始值，不做滤波处理
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);
double t1 = cv::getTickCount();
cout << "time=" << (t1 - t0) / cv::getTickFrequency() << endl;
imwrite("data/test/srcImg.jpg", srcImg);
imwrite("data/test/myFilter.jpg", dstImg);
}
#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);
}
参考：https://www.cnblogs.com/riddick/p/7989871.html
扫一扫
1231
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
