更新式中值滤波算法

注：目前开通个人网站朝思录，之后的博文将在上面更新，CSDN博客会滞后一点

---

前言

中值滤波算法是一种非线性的滤波算法，其中心思想是采用模板内所有像素的排序中值作为目标像素的值，进行滤波。一般情况下，若模板大小为$m\times m$，可有效滤除面积小于$m^2/2$的脉冲像素团。比如$3\times3$中值滤波模板可有效滤除面积为1的椒盐噪声。
最近在写中值滤波算法时，因为一些bug偶然发现了一些有趣的现象。中值滤波算法的大致算法流程为：

% 算法1
% 读取图像
im = imread('某个地址');
% 中值滤波
im_backup = im;     %备份图像
for i = size(im_backup,1)
    for j = size(im_backup,2)
        value = GetMedian(im_backup,i,j,3);     %获取3*3模板的中值
        im(i,j) = value;
    end
end

注：GetMedian为自定义的获取模板中值的函数。
其中比较关键的部分是备份图像，得到im和im_backup，并采用im_backup逐像素得出中值修改im。备份图像的目的是防止已修改的像素对下一像素的模板中值产生影响。（在线性滤波中这种备份很重要，否则滤波就不再具有线性性）。
但是，刚开始写程序时，笔者忘记备份图像，此时笔者算法流程如下所示：

% 算法2
% 读取图像
im = imread('某个地址');
% 中值滤波
for i = size(im,1)
    for j = size(im,2)
        value = GetMedian(im,i,j,3);        %获取3*3模板的中值
        im(i,j) = value;
    end
end

但是采用该“错误”的算法取得的结果出乎意料的好，椒盐噪声滤除非常干净，甚至比原来的中值滤波算法好。因为算法2中，除了第一次迭代，之后所有迭代中模板里总包含被更新的像素。为了方便区分，故称算法2为更新式中值滤波算法，算法1为传统中值滤波算法。

对比

1、椒盐噪声滤除效果的对比

原始图像采用冈萨雷斯第三版中受椒盐噪声污染的PCB图。

采用的模板大小为$3\times3$，对比结果如表所示：

$3\times3$传统中值滤波算法	$3\times3$更新式中值滤波算法

可以发现在同样的滤波条件下，更新式中值滤波算法在牺牲轻微锐度（如QFP下的总线边缘略不清晰）的基础上滤除了更多噪声。而牺牲的锐度可以通过锐化来补偿。
当然，在模板大小为$5\times5$时，传统中值滤波算法也会滤除掉大部分噪声，但此时图像引入了非常多的blurring，导致图像变得更加模糊。对比如下图所示。

$5\times5$传统中值滤波	$3\times3$更新中值滤波

所以在某些应用中，可能采用更新式中值滤波算法效果较好。除了椒盐噪声的去除效果更好之外，其不需要备份图像更节省了存储空间。

2、彩色图像滤波效果对比

待滤波图像如下图所示。该图一度作为笔者的壁纸使用。出处不详，往知情者告知。

这里都采用$7\times7$模板进行滤波。滤波结果图如下所示。

传统中值滤波算法

更新式中值滤波算法

可以看到更新式中值滤波算法滤除了更多的高频信息，如图片左侧中部的梯田：传统中值滤波算法导致了不自然的高频横纹；而更新式滤波算法并没有相应的纹路。而且更新式中值滤波算法得到的图像色块偏浓重，略有野兽派的画风；传统中值滤波算法在近景出还有一定的写实特点。

附代码

脚本

1、Median_Filter_Script.m

clc
clear
close all

%% 定义
% 滤波器模板大小
s = 5;

% 取值位置
% 注（中值滤波为0.5，最大值滤波为1，最小值滤波为0）
p = 0.5;

% 更新式滤波
% 1：非更新式滤波，0：更新式滤波
method = 1;
%% 图像输入
im = imread('circuit_board.tif');
im = im(:,:,1);

%% 图像扩展
% 根据滤波器模板大小定义扩展图像
% Define augmented figure according to the filter mask scale
[H,W] = size(im);
[im_augmented, index_x, index_y] = Image_Augment(im, s);

%% 图像滤波
% 注意，如果采用不更新方式进行滤波，出来的效果貌似更好
im_filtered_augmented = Median_Filter(im_augmented,s,p,method);
im_filtered = im_filtered_augmented(index_x,index_y);

%% 图像显示
figure
subplot(1,2,1
)
imshow(im);
subplot(1,2,2)
imshow(uint8(im_filtered));

2、Color_Image_Median_Filter_Script.m

clc
clear
close all
%% 定义
% 滤波器模板大小
s = 7;

% 取值位置
% 注（中值滤波为0.5，最大值滤波为1，最小值滤波为0）
p = 0.5;

% 更新式滤波
% 1：非更新式滤波，0：更新式滤波
method = 0;
%% 图像输入
im = imread('test.jpg');

%% 图像增广
[im_R_augmented,index_R_x,index_R_y] = Image_Augment(im(:,:,1),s);
[im_G_augmented,index_G_x,index_G_y] = Image_Augment(im(:,:,2),s);
[im_B_augmented,index_B_x,index_B_y] = Image_Augment(im(:,:,3),s);

%% 滤波
im_filtered_R_augmented = Median_Filter(im_R_augmented,s,p,method);
im_filtered_G_augmented = Median_Filter(im_G_augmented,s,p,method);
im_filtered_B_augmented = Median_Filter(im_B_augmented,s,p,method);

im_filtered = cat(3,...
    im_filtered_R_augmented(index_R_x,index_R_y),...
    im_filtered_G_augmented(index_G_x,index_G_y),...
    im_filtered_B_augmented(index_B_x,index_B_y));

%% 显示图像
imshow(uint8(im_filtered));

函数

3、Image_Augment.m

function [im_augmented,index_h,index_w] = Image_Augment(im,s)
% Image_Augment - 增广图像，以便于图像模板匹配
% 输入参数
% im - 待增广图像
% s - 模板边长
% 输出参数
% im_augmented - 增广后图像
% index_h,index_w - 图像在高和宽的嵌入位置索引

[H,W] = size(im);
index_h = s:(s+H-1);
index_w = s:(s+W-1);
im_augmented = zeros(H + 2*(s-1),...
    W + 2*(s-1));
im_augmented(index_h,index_w) = double(im);

im_augmented(index_h,1:(s-1)) = repmat(...
    im_augmented(index_h,s),1,s-1);
im_augmented(index_h,(s+W):(W+2*(s-1))) = repmat(...
    im_augmented(index_h,s+W-1),1,s-1);
im_augmented(1:(s-1),index_w) = repmat(...
    im_augmented(s,index_w),s-1,1);
im_augmented((s+H):(H+2*(s-1)),index_w) = repmat(...
    im_augmented(s+H-1,index_w),s-1,1);

for i = 1:(s-1)
    for j = 1:(s-1)

        if i > j
            im_augmented(i,j) = im_augmented(s,j);
            im_augmented(H+2*(s-1)-i+1,j) = ...
                im_augmented(s+H-1,j);
            im_augmented(i,W+2*(s-1)-j+1) = ...
                im_augmented(s,W+2*(s-1)-j+1);
            im_augmented(H+2*(s-1)-i+1,W+2*(s-1)-j+1) = ...
                im_augmented(s+H-1,W+2*(s-1)-j+1);
        elseif i < j
            im_augmented(i,j) = im_augmented(i,s);
            im_augmented(H+2*(s-1)-i+1,j) = ...
                im_augmented(H+2*(s-1)-i+1,s);
            im_augmented(i,W+2*(s-1)-j+1) = ...
                im_augmented(i,s+W-1);
            im_augmented(H+2*(s-1)-i+1,W+2*(s-1)-j+1) = ...
                im_augmented(H+2*(s-1)-i+1,s+W-1);
        else
            im_augmented(i,j) = im_augmented(s,s);
            im_augmented(H+2*(s-1)-i+1,j) = ...
                im_augmented(s+H-1,s);
            im_augmented(i,W+2*(s-1)-j+1) = ...
                im_augmented(s,s+W-1);
            im_augmented(H+2*(s-1)-i+1,W+2*(s-1)-j+1) = ...
                im_augmented(s+H-1,s+W-1);
        end
    end        
end
end

4、Median_Filter.m

function im_filtered = Median_Filter(im,s,ratio,method)
% Median_Filter - 图像中值滤波
% 输入参数
% im - 待滤波图像
% s - 滤波模板大小
% ratio - 滤波阈值
% method - 滤波方法
%   1：非更新式滤波，0：更新式滤波
%   注：非更新式滤波最好只采用0.5滤波阈值，采用其他阈值会出现不理想的效果
% 输出参数
% im_filtered - 滤波后图像

%% 基本参数求解
[H,W] = size(im);

%% 图像增广
[im_augmented,index_x,index_y] = Image_Augment(im,s);
im_augmented_temp = im_augmented;
for i = (s-1):(H+2*(s-1)-(s-2))
    for j = (s-1):(W+2*(s-1)-(s-2))
        if method
            temp = im_augmented_temp(...
                (i-(s-2)):(i+(s-2)),(j-(s-2)):(j+(s-2)));
        else
            temp = im_augmented(...
                (i-(s-2)):(i+(s-2)),(j-(s-2)):(j+(s-2)));
        end
        temp = sort(temp(:));
        index = ceil(length(temp)*ratio);
        if index == 0
            index = 1;
        end
        im_augmented(i,j) = temp(index);
    end
end
im_filtered = im_augmented(index_x,index_y);
end