上一篇分享了一些均值滤波相关的算法，均值滤波作为一种线性滤波器，在滤除噪声的同时也会导致边缘模糊问题。而且均值滤波对高斯噪声的效果很好，但是对于椒盐噪声的效果就很一般。但是中值滤波作为一种顺序滤波器，对于椒盐噪声的效果很好，而且保边能力很强，所以这一篇主要讨论一下中值相关的算法。

中值滤波

算法原理

中值滤波很好理解，均值滤波就是在一个小窗口中求均值来取代当前像素值，而中值滤波就是通过求小窗口中的中位值来取代当前位置的方式来滤波。

如图绿色窗口就是当前的滤波窗口，在一个3X3的邻域窗口中进行滤波。那么中值滤波做的就是：

1. 对这个邻域中的像素值进行排序：[45, 50, 52, 60, 75, 80, 90, 200, 255];
2. 从排序后的数据中找出中位值：75；
3. 用中位值取代当前位置的像素值，所以就可以得到右侧的滤波后的数据。

代码实现

%% --------------------------------
%% author:wtzhu
%% email:wtzhu_13@163.com
%% date: 20210306
%% fuction: 中值滤波
%% --------------------------------
clear;
clc;
close all;
img = imread('./images/test_pattern_blurring_orig.tif');
[m, n] = size(img);
figure;subplot(221);imshow(img);title('original image');

%% 运算的时候需要对边缘进行扩展
% 需要特殊处理四周最外圈的行和列，本算法中将其向外扩展一圈，用最外圈的值填充
headRowMat = img(1,:);%取f的第1行
tailRowMat = img(m,:);%取f的第m行
% 行扩展后，列扩展时需要注意四个角需要单独扩展进去，不然就成了十字架形的
headColumnMat = [img(1,1), img(:,1)', img(m,1)];
tailColumnMat = [img(1,n), img(:,n)', img(m,n)];
expandImage = [headRowMat; img; tailRowMat];
expandImage = [headColumnMat; expandImage'; tailColumnMat];
expandImage = uint8(expandImage');
subplot(222);imshow(expandImage);title('expand image');

newImg = zeros(m, n);
for i =2: m+1
    for j =2: n+1
       imgRoi = [expandImage(i-1, j-1) expandImage(i-1, j) expandImage(i-1, j+1) ...
                 expandImage(i  , j-1) expandImage(i  , j) expandImage(i  , j+1) ...
                 expandImage(i+1, j-1) expandImage(i+1, j) expandImage(i+1, j+1)];
       orderedList = sort(imgRoi);
       sizeRoi = size(imgRoi);
       newImg(i-1, j-1) = orderedList((sizeRoi(2)+1)/2);
    end
end
newImg = uint8(newImg);
subplot(223);imshow(newImg);title('new image');
subplot(224);imshow(uint8(newImg-img));title('newImg-img');

从滤波效果图中可以看出去噪能力还可以，百块中的噪点去除了很多，而且边缘信息保留得也很好。

多级中值滤波

算法原理

简单的中值滤波器滤波效果有限，于是就有人提出了将多个中值滤波进行多级级联实现更好的滤波效果。

如图就是一种多级级联的方式，先在窗口中定义一个’+'和’X’形的窗口，然后分别求出这两个窗口的中位值，然后结合当前窗口的中心点就有3个候选值，再从这三个值中求出一个中位值作为滤波后的结果。

这种方式也可以直接应用到RAW图中做BNR，需要修改的就是窗口设置为5X5，然后在做滤波的时候需要区分G和RB通道，因为这个前面已经讲过，RAW图中的RGB分布是不均匀的，G占50%，R和B各占25%。

左侧就是针对G通道的滤波器，右侧是R和B通道的滤波器，都是定义了一个’+'和’X’形的窗口，不同的只是取的点的位置不同。

代码实现

%% --------------------------------
%% author:wtzhu
%% email:wtzhu_13@163.com
%% date: 20211023
%% fuction: multistage median filters
%% --------------------------------

close all;
clear all;
clc
img = imread('./images/lena.bmp');
I = double(img);
% I = double(imresize(img, [64, 64]));
figure();
imshow(uint8(I));
title('org file');

I_noise = I + 10 * randn(size(I));
figure();
imshow(uint8(I_noise));
title('noise file');
[m,n] = size(I_noise);

DenoisedImg = zeros(m,n);
PaddedImg = padarray(I,[1, 1],'symmetric','both');

tic
for i = 1: m
    for j = 1: n
        roi = PaddedImg(i:i+2, j:j+2);
        % first stage
        median_HV = median([roi(1,2), roi(2,1), roi(2,2), roi(2,3), roi(3,2)]);
        median_diag = median([roi(1,1), roi(1,3), roi(2,2), roi(3,1), roi(3,3)]);
        % second stage
        DenoisedImg(i, j) = median([median_HV, roi(2,2), median_diag]);
    end
end
figure();
imshow(uint8(DenoisedImg));
title('denoise file');
toc

b = medfilt2(I_noise,[3,3]);
figure();
imshow(uint8(b));
title('median filter of matlab denoise file');

多级中值混合滤波

算法原理

前面介绍过中值滤波和均值滤波各自有各自的优点和缺点，所以就可以考虑将两者结合起来，相互弥补实现更好的滤波效果，于是就有人提出了这种多级混合滤波的方式。

算法流程：

1. 求出竖直方向相邻三个点的均值和水平方向相邻三个点的均值，再结合当前点，用这三个点再求一个中位值；
2. 求出45°和135°方向上的均值，然后结合当前点求出一个中位值；
3. 两个中位值结合当前点组成新的数组，最后求一个中位值作为当前点的值完成滤波。

这样就巧妙的将均值滤波和中位值滤波结合了。然后如果应用在RAW图上，只需要对滤波器稍作改善即可。

代码实现

%% --------------------------------
%% author:wtzhu
%% email:wtzhu_13@163.com
%% date: 20211023
%% fuction: multistage median hybird filters
%% --------------------------------

close all;
clear all;
clc
img = imread('./images/lena.bmp');
I = double(img);
figure();
imshow(uint8(I));
title('org file');

I_noise = I + 10 * randn(size(I));
figure();
imshow(uint8(I_noise));
title('noise file');
[m,n] = size(I_noise);

DenoisedImg = zeros(m,n);
PaddedImg = padarray(I,[1, 1],'symmetric','both');

tic
for i = 1: m
    for j = 1: n
        roi = PaddedImg(i:i+2, j:j+2);
        % first stage: average and median
        mean_V = mean(roi(:,2));
        mean_H = mean(roi(2,:));
        median_HV = median([mean_V, roi(2, 2)], mean_H);
        
        mean45 = mean([roi(1, 3), roi(2, 2), roi(3, 1)]);
        mean135 = mean([roi(1, 1), roi(2, 2), roi(3, 3)]);
        median_diag = median([mean45, roi(2, 2)], mean135);
        
        % second stage
        DenoisedImg(i, j) = median([median_HV, roi(2,2), median_diag]);
    end
end
figure();
imshow(uint8(DenoisedImg));
title('denoise file');
toc

b = medfilt2(I_noise,[3,3]);
figure();
imshow(uint8(b));
title('median filter of matlab denoise file');

多级中值有理混合滤波

算法原理

WMF加权中值滤波

加权中值滤波也很好理解，和加权均值滤波差不多，就是在原始数据的基础上给每个点分别赋予一个权重，然后在加权后的数据中取出中位值作为滤波后的值。

算法流程

1. 求出’+'形和’X’形的窗口的中位值；

2. 对’+'形窗口再利用CWMF求出一个值，CWMF是WMF的一种特殊情况，就是只对中心点进行加权；

3. 对以上求出的三个参数用一下公式计算出一个新的值作为滤波后的值
   $$y(m,n)=\varphi 2(m,n)+\frac{\varphi 1(m,n)-2\ast \varphi 2(m,n)+\varphi 3(m,n)}{h+k(\varphi 1(m,n)-\varphi 3(m,n))}$$

一样的，稍作改动该算法就可以用于raw格式图像

代码实现

%% --------------------------------
%% author:wtzhu
%% email:wtzhu_13@163.com
%% date: 20211023
%% fuction: median rational hybird file filter
%% references: MEDIAN-RATIONAL HYBRID FILTERS
%% --------------------------------

close all;
clear all;
clc
img = imread('./images/lena.bmp');
I = double(img);
figure();
imshow(uint8(I));
title('org file');

I_noise = I + 10 * randn(size(I));
figure();
imshow(uint8(I_noise));
title('noise file');
[m,n] = size(I_noise);

DenoisedImg = zeros(m,n);
PaddedImg = padarray(I,[1, 1],'symmetric','both');

h = 2;
k = 0.01;

tic
for i = 1: m
    for j = 1: n
        roi = PaddedImg(i:i+2, j:j+2);
        median_HV = median([roi(1,2), roi(2,1), roi(2,2), roi(2,3), roi(3,2)]);
        median_diag = median([roi(1,1), roi(1,3), roi(2,2), roi(3,1), roi(3,3)]);
        CWMF = median([roi(1,2), roi(2,1), roi(2,2)*3, roi(2,3), roi(3,2)]);
        
        DenoisedImg(i, j) = CWMF + (median_HV + median_diag - 2 * CWMF) / (h + k * (median_HV - median_diag));
    end
end
toc
figure();
imshow(uint8(DenoisedImg));
title('denoise file');
b = medfilt2(I_noise,[3,3]);
figure();
imshow(uint8(b));
title('median filter of matlab denoise file');

相关链接

• zhihu： ISP图像处理 - 知乎 (zhihu.com)
• CSDN：ISP图像处理_wtzhu_13的博客-CSDN博客
• Bilibili：食鱼者的个人空间_哔哩哔哩_Bilibili
• Gitee：ISPAlgorithmStudy: ISP算法学习汇总，主要是论文总结 (gitee.com)

一、自适应中值滤波（代码）

function aOut = stPlus(aIn,f1,f2,st,stMax)
%这仅仅是 AdaMedFilter (自适应中值滤波函数) 的子函数
%这仅仅是一个附属函数/功能片段（有特定的功能，输入若干边界条件，输出一个具体的数值）
%输入aIn，需要滤波的图像，循环到的行数f1和列数f2，掩模版大小st，和掩模版最大值stMax

[ra,ca] = size(aIn);
%制作出掩模版
sxy = aIn(f1-round(st/2)+1:f1+round(st/2)-1,f2-round(st/2)+1:f2+round(st/2)-1);
%将掩模版从小到大排序并输出为一列
zSort = sort(sxy(:));
%基于排序后的zSort找到最大值，最小值，中值，以及循环位置的值
med = floor(st*st/2)+1;
zMax = zSort(st*st);
zMin = zSort(1);
zMed = zSort(med);
zxy = sxy(floor(st/2)+1,floor(st/2)+1);
%判断中值是否在最大值和最小值之间
delta1 = zMed-zMin;
delta2 = zMax-zMed;
if (delta1 * delta2) > 0
    %判断循环位置的值是否在最大值和最小值之间
    delta3 = zxy - zMin;
    delta4 = zMax - zxy;
    if (delta3 * delta4) > 0
        %中值 在 最大最小值之间，且循环位置的值 在 最大最小值之间
        aOut = zxy;
    else
         %中值 在 最大最小值之间，但循环位置的值 不在 最大最小值之间
        aOut = zMed;
    end
end
if (delta1 * delta2) <= 0 
    %中值 不在 最大最小值之间
    if st < stMax
        %掩模版尺寸 不是 最大
        if f1 >= floor((st+2)/2)+1 && f2 >= floor((st+2)/2)+1 && f1 <= ra-floor((st+2)/2) && f2 <= ca-floor((st+2)/2)
        %该循环位置是否还可以进行st+2的循环
            aOut = Inf;%可以
        else
            aOut = zMed;%不可以，直接输出
        end
    else
        %掩模版尺寸 是 最大，直接输出
        aOut = zMed;
    end
end
end

function aOut = AdaMedFilter(aIn,static,stMax)
%这仅仅是一个 自适应 中值滤波函数
%输入aIn图像，输出aOut结果
%static是人为定制的静态变量，意义是sxy掩模版的大小，默认为3*3，stMax就是掩模版最大的大小，默认为7*7

%输入参数的默认及处理
if nargin == 1
    static = 3;
    stMax = 7;
elseif nargin == 2
    stMax = 7;
end

%定义输入及输出图像
%aIn是最先输入图像
%atemp是对输入图像的初步处理（比如归一化，双精度化，归255化等），可以在后面调用时灵活更改
%aOut是最终输出图像
atemp = aIn;
[ra,ca] = size(aIn);
aOut = zeros(ra,ca);

%定义st = 3，其中static为人为定义的静态变量，所以
st = static;

%循环，每次循环把循环位置的值给进去，输出一个唯一的有效值赋值给aOut在当前位置的值
for f1 = round(st/2):ra-round(st/2)+1
    for f2 = round(st/2):ca-round(st/2)+1
        %st = static在每一个循环的循环开始和循环结束都有，是为了绝对保证st在每次循环开始和末尾都是绝对为3
        st = static;
        %输入循环位置的值，输出处理后的值，当然，有可能st = 3的中值滤波不够，那么就会输出temp = Inf
       temp = stPlus(atemp,f1,f2,st,stMax);
       %我已经在函数stPlus中确定当temp = Inf时，一定能够被st+2的中值滤波处理，所以我就循环让st增大取得有效的temp
        while(temp == Inf)
           st = st+2;
           temp = stPlus(atemp,f1,f2,st,stMax);
        end
        %每一个循环最重要给我一个唯一的正确值temp，我把它赋值给循环所在位置的aOut！
        aOut(f1,f2) = temp;
        %st = static在每一个循环的循环开始和循环结束都有，是为了绝对保证st在每次循环开始和末尾都是绝对为3
        st = static;
    end
end
end

clear
close all
aIn = imread("GaosiFilter.jpg");
static = [3,7,11];
stMax = 11;

%取不同的st值（掩模版大小）进行滤波
a1 = AdaMedFilter(aIn,static(1),stMax);
a2 = AdaMedFilter(aIn,static(2),stMax);
a3 = AdaMedFilter(aIn,static(3),stMax);

%对滤波后图片a1进行均值滤波以消除图片轮廓的“毛刺”
a1Ave = AverageFilter(a1,5);
figure,imshow(a1Ave,[])

%反复进行滤波3次
a12 = AdaMedFilter(a1,static(1),stMax);
a13 = AdaMedFilter(a12,static(1),stMax);
a14 = AdaMedFilter(a13,static(1),stMax);
a15 = AdaMedFilter(a14,static(1),stMax);

figure,imshow(aIn,[])
figure,imshow(a1,[])
figure,imshow(a2,[])
figure,imshow(a3,[])

figure,imshow(a1,[])
figure,imshow(a12,[])
figure,imshow(a13,[])
figure,imshow(a14,[])
figure,imshow(a15,[])

%对比1次滤波和5次滤波的剖线差别
figure,hold on;
plot(a1(220,:),'r')
plot(a15(220,:),'b')

二、自适应中值滤波（图片）

1、输入的噪声图像

2、当st = 3，stmax = 11时输出的滤波后图像

3、分析

中值滤波对胡椒噪声和盐噪声滤波效果特别好，但是自适应保边中值滤波会让滤波后图片的边缘出现“毛刺”，这个时候我们可以采用均值滤波，以图片稍微变模糊一些的代价将边缘的毛刺消除，如下图：

三、st的不同对滤波效果的影响

1、st = 3

2、st = 7

3、st = 11

4、结论

​ 随着st（掩模版大小）的增大，滤波效果并没有变得更好，甚至会让原本清晰的图片边缘变得更模糊。

四、探究对图像重复进行滤波多次对滤波效果的影响

1、滤波一次图片

2、滤波五次图片

3、对比一次滤波的剖线和五次滤波的剖线

说明：

蓝色线条是“5次滤波”图片的剖线；红色线条是“1次滤波”图片的剖线

4、结论

​ 从图片上看，一次滤波和五次滤波效果差别不大，但是从剖线上看，很明显的感知到，对于图片平滑部位，一次滤波和五次滤波效果相差无几，但是对于图像急剧变化的部分，很明显的是五次滤波后的图像线条更加的平滑，变化更加均匀！