直方图实现快速中值滤波

中值滤波能够有效去除图像中的异常点，具有去除图像噪声的作用。传统中值滤波的算法一般都是在图像中建立窗口，然后对窗口内的所有像素值进行排序，选择排序后的中间值作为窗口中心像素滤波后的值。由于这个做法在每个像素点处都要建立窗口并排序，非常耗时，尤其是有大量的冗余计算。如下图：

黄色区域+中间粉色区域是第一个像素为中心建立的滤波窗口，粉色区域+右边蓝色区域为同一行第二个像素为中心建立的滤波窗口。传统做法对每一个窗口内所有像素排序，那么中间粉色区域的像素就会被重复的做排序运算，存在大量冗余计算。如果窗口移动一个像素的时候只用考虑去除左边一列内（黄色区域）的像素，并且加上右边一列（蓝色区域）的像素，运算会大大简化。这样的操作可以使用直方图来实现。

一、直方图实现快速中值滤波算法流程：

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 }
View Code

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 }
View Code
posted @ 2017-12-05 22:42  一度逍遥  阅读(7663)  评论(4)  编辑  收藏  举报