直方图实现快速中值滤波

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

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

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

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:

#include <opencv2opencv.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);
}

Code02:

#include <opencv2opencv.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);
}