图像处理——（源）中值滤波（MediFilter）函数编程实现

https://blog.csdn.net/tengfei461807914/article/details/83626123

中值滤波是一种非线性滤波，在处理脉冲噪声以及椒盐噪声时效果极佳，能够有效的保护好图像的边缘信息。

中值滤波的处理思路很简单，取卷积核当中所覆盖像素中的中值作为锚点的像素值即可。

如果按照遍历所有像素，再对卷积核中的像素排序取中值，那么时间复杂度会很高，需要对中值滤波进行改进。

中值滤波的改进实际上很是很好想的，无非就是一个滑动窗口取中值的问题，每次向右滑动的过程中等于在窗口中新添加添加一列窗口像素，同时减去一列窗口像素，考虑维护这个窗口中的像素信息变化即可。

这里面使用huang算法，该思路是这个人提出来的，算法思路如下：

在计算中值的办法中，不使用排序，而是使用像素直方图，也就是记录像素值的哈希。首先设定阈值threshold，这个threshold就是窗口的中心位置，即ksize×ksize/2+1，kisze为窗口尺寸。

每次在计算中值的过程中，从小到大累加像素直方图的值，如果该值大于等于，此时对应的像素值就是中值了。

例如ksize=3的窗口如下：
⎡⎣⎢122235154⎤⎦⎥(3) left[egin{matrix}1 & 2 & 1 \2 & 3 & 5 \2 & 5 & 4end{matrix} ight] ag{3}
⎣
⎡
​

1
2
2
​

2
3
5
​

1
5
4
​

⎦
⎤
​
(3)

对该窗口中的值计算像素直方图如下，threshold=3×3/2+1=5
1:2（表示像素值为1的有2个）
2:3
3:1
4:1
5:2

因为2+3≥5，所以中值为2

每次滑动窗口的过程中，如果窗口在第一列，那么直接计算直方图。否则向右移动，在直方图中减去左侧离开窗口中像素，在右侧添加进入窗口中的像素。

此外，也可以让窗口按照蛇形来移动，这样也会避免每次窗口在第一列时需要重新计算的问题。

void AddSultPapperNoise(const Mat &src, Mat &dst,int num)//添加椒盐噪声
{
    dst = src.clone();
    uchar *pd=dst.data;
    int row, col, cha;
    srand((unsigned)time(NULL));
    while (num--)
    {
        row = rand() % dst.rows;
        col = rand() % dst.cols;
        cha = rand() % dst.channels();
        pd[(row*dst.cols + col)*dst.channels() + cha] = 0;
    }
}


int GetMediValue(const int histogram[], int thresh)//计算中值
{
    int sum = 0;
    for (int i = 0; i < (1 << 16); i++)
    {
        sum += histogram[i];
        if (sum >= thresh)
            return i;
    }
    return (1 << 16);
}

void MyFastMediFilter(const Mat &src, Mat &dst, int ksize)
{
    CV_Assert(ksize % 2 == 1);

    Mat tmp;
    int len = ksize / 2;
    tmp.create(Size(src.cols + len, src.rows + len), src.type());//添加边框
    dst.create(Size(src.cols, src.rows), src.type());


    int channel = src.channels();
    uchar *ps = src.data;
    uchar *pt = tmp.data;
    for (int row = 0; row < tmp.rows; row++)//添加边框的过程
    {
        for (int col = 0; col < tmp.cols; col++)
        {
            for (int c = 0; c < channel; c++)
            {
                if (row >= len && row < tmp.rows - len && col >= len && col < tmp.cols - len)
                    pt[(tmp.cols * row + col) * channel + c] = ps[(src.cols * (row - len) + col - len) * channel + c];
                else
                    pt[(tmp.cols * row + col) * channel + c] = 0;
            }
        }
    }
    int Hist[(1 << 16)] = { 0 };
    uchar *pd = dst.data;
    ushort val = 0;
    pt = tmp.data;
    for (int c = 0; c < channel; c++)//每个通道单独计算
    {
        for (int row = len; row < tmp.rows - len; row++)
        {
            for (int col = len; col < tmp.cols - len; col++)
            {
                
                if (col == len)
                {
                    memset(Hist, 0, sizeof(Hist));
                    for (int x = -len; x <= len; x++)
                    {
                        for (int y = -len; y <= len; y++)
                        {
                            val = pt[((row + x) * tmp.cols + col + y) * channel + c];
                            Hist[val]++;
                        }
                    }
                }
                else
                {
                    int L = col - len - 1;
                    int R = col + len;
                    for (int y = -len; y <= len; y++)
                    {
                        int leftInd = ((row + y) * tmp.cols + L) * channel + c;
                        int rightInd = ((row + y) * tmp.cols + R) * channel + c;
                        Hist[pt[leftInd]]--;
                        Hist[pt[rightInd]]++;
                    }
                }
                val = GetMediValue(Hist, ksize*ksize / 2 + 1);
                pd[(dst.cols * (row - len) + col - len) * channel + c] = val;
                
            }
        }
    }
}