基于MATLAB的数字图像处理（二）——空间域滤波篇（不使用内置函数）

目录

概述

空间域平滑——低通滤波

空间域锐化——高通滤波

---

概述

图像的滤波分为空间域和频率域两方面，在此首先说明空间域滤波。空间域滤波又分为高通滤波和低通滤波两种，下面分别阐述。

空间域平滑——低通滤波

低通滤波可以理解为降噪的过程，常用的方法有均值滤波、高斯滤波、中值滤波、K个近邻的平滑、梯度倒数加权平滑等。

均值滤波，顾名思义即是取范围内各点灰度的平均值作为中间像元的灰度值来进行滤波，3*3模板的演示代码如下：

%%
% 均值滤波
clear;
clc;
%%
% 主程序入口
Main();
%%
% 主函数
function Main()
    Image = imread('lena.tif');
    NoiseImage = imnoise(Image,'salt & pepper',0.02);
    OutImage = Method_GaussFilter(NoiseImage);
    figure(1);
    imshow(NoiseImage);
    title('椒盐噪声图像');
    figure(2);
    imshow(OutImage);
    title('均值滤波图像');
end
%%
% 均值滤波函数
function OutImage = Method_GaussFilter(NoiseImage)
    OutImage = double(NoiseImage);
    [Height,Width] = size(NoiseImage);
    Temp = ones(3)/9;
    % Temp为3*3均值滤波模板
    for i = 2:Height-1
        for j = 2:Width-1
            ImageTemp = NoiseImage(i-1:i+1,j-1:j+1);
            OutImage(i,j) = ConvImage(Temp,ImageTemp);
        end
    end
    OutImage = uint8(OutImage);
end
%%
% 卷积函数
function OutTemp = ConvImage(Temp,ImageTemp)
    OutTemp = sum(sum(Temp.*double(ImageTemp)));
end

高斯低通滤波可以理解为在均值滤波的基础上对像元进行了加权，只不过权满足一定的函数规律，例如中央像元权重最大等，示例代码如下：

%%
% 高斯滤波
clear;
clc;
%%
% 主程序入口
Main();
%%
% 主函数
function Main()
    Image = imread('lena.tif');
    NoiseImage = imnoise(Image,'salt & pepper',0.02);
    OutImage = Method_GaussFilter(NoiseImage);
    figure(1);
    imshow(NoiseImage);
    title('椒盐噪声图像');
    figure(2);
    imshow(OutImage);
    title('高斯滤波图像');
end
%%
% 高斯滤波函数
function OutImage = Method_GaussFilter(NoiseImage)
    OutImage = double(NoiseImage);
    [Height,Width] = size(NoiseImage);
    Temp = zeros(3);
    sigma = 0.8;
    for i=1:3
        for j=1:3
            Temp(i,j)=exp(-((i-1-1)^2+(j-1-1)^2)/2/sigma^2)/(2*pi*sigma^2);
            % Temp为高斯模板
        end
    end
    for i = 2:Height-1
        for j = 2:Width-1
            ImageTemp = NoiseImage(i-1:i+1,j-1:j+1);
            OutImage(i,j) = ConvImage(Temp,ImageTemp);
        end
    end
    OutImage = uint8(OutImage);
end
%%
% 卷积函数
function OutTemp = ConvImage(Temp,ImageTemp)
    OutTemp = sum(sum(Temp.*double(ImageTemp)));
end

下面介绍最好用的一种滤波器——中值滤波器，通过代码运行即可看到，均值滤波的效果明显优于以上两种，原因则是它的滤波方式并不是卷积，而是对中央像元作了中位数替代，从而并没有以牺牲其他非噪声像元的真实值为代价。示例代码如下：

%%
% 中值滤波
clear;
clc;
%%
% 主程序入口
Main();
%%
% 主函数
function Main()
    Image = imread('lena.tif');
    NoiseImage = imnoise(Image,'salt & pepper',0.02);
    OutImage = Method_MedFilter(NoiseImage);
    figure(1);
    imshow(NoiseImage);
    title('椒盐噪声图像');
    figure(2);
    imshow(OutImage);
    title('中值滤波图像');
end
%%
% 中值滤波函数
function OutImage = Method_MedFilter(NoiseImage)
    [Height,Width] = size(NoiseImage);
    OutImage = zeros(Height,Width);
    for i = 2:Height-1
        for j = 2:Width-1
            OutImage(i,j) = SortImage(NoiseImage(i-1:i+1,j-1:j+1));
        end
    end
    OutImage = uint8(OutImage);
end
%%
% 排序筛选函数
function OutTemp = SortImage(Temp)
    Temp = sort(Temp);
    OutTemp = Temp(5);
end

下面则是K个近邻平滑，相对于梯度倒数加权较容易理解，它首先选取中央像元一定邻域内的像元逐个求差，得到的差的绝对值最小的n个对应的像元再求平均。对于3*3邻域，n取6为最佳。示例代码如下：

%灰度最接近的K个邻点平滑，MATLAB实现
%%
% 主程序
clear;
clc;

SizeTemp = 3;
Dich = (SizeTemp-1)/2;
Image = imread('lena.tif');  
figure(1);
imshow(Image);
ImageSP = imnoise(Image,'salt & pepper',0.02); 
figure(2);
imshow(ImageSP);
[Height,Width] = size(ImageSP);
ImageSP = double(ImageSP);
MedImage = ImageSP;
for k = 1:6
    for i = Dich+1:Height-Dich
        for j = Dich+1:Width-Dich
            m = ImageSP(i-Dich:i+Dich,j-Dich:j+Dich) - ImageSP(i,j);
            m(Dich+1,Dich+1) = 256;
            Temp = absmin(m(:));
            MedImage(i,j) = ImageSP(i,j) + sum(Temp(1:k))/k;
        end
    end
end
figure(3);
imshow(uint8(MedImage));
title(strcat('k = ',num2str(k)))

function out=absmin(m)
    for i=1:size(m)
        for j=1:size(m)-i
            if(abs(m(j))>abs(m(j+1)))
                temp=m(j);
                m(j)=m(j+1);
                m(j+1)=temp;
            end
        end
    end
out=m;
end

需要注意的是此段代码运行时间较长，因此推测存在可优化空间，此处不再赘述，下面讨论梯度倒数加权平滑，因为相对较难理解，此处直接上代码：

%梯度倒数加权平滑
%%
% 主程序
Image = imread('lena.tif');
ImageSP = imnoise(Image,'salt & pepper',0.02);
OutImage = Grs(ImageSP);
imshow(OutImage);

function OutImage = Grs(Image)
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
    %%
    % 基本参数及预处理
    [Height,Width]=size(Image);
    Image = double(Image);
    SizeTemp = 3;
    Temp = zeros(SizeTemp);
    %%
    % 构造模板
    for i = 2:Height-1
        for j = 2:Width-1
            Temp=zeros(3);
            for m=-1:1
               for n=-1:1
                  if(Image(i+m,j+n)-Image(i,j)==0)
                      Temp(m+2,n+2)=0;
                  else
                      Temp(m+2,n+2)=1/abs(Image(i+m,j+n)-Image(i,j));
                  end
               end
           end
        end
    end
    %%
    % 卷积
    OutImage = uint8(double(conv_image(Image,Temp)));
end

空间域锐化——高通滤波

高通滤波与边缘提取有一定相似之处。常见的有Robert算子、Sobel算子、Laplace算子等，下面主要针对Laplace算子进行讨论。

通过代码可以看出锐化算子与边缘提取算子的不同之处在于锐化算子是针对边缘进一步加深了梯度，而边缘检测算子仅仅是提取出边缘

锐化：

%%
% laplace算子
clear;
clc;
%%
% 主程序入口
Main();
%%
% 主函数
function Main()
    Image = imread('lena.tif');
    OutImage = uint8(Method_Laplace(Image));
    figure('name','Laplace锐化');
    imshow(OutImage);
    title('Laplace锐化图像');
end
%%
% laplace算子函数
function OutImage = Method_Laplace(Image)
    w1 = [0 -1 0;-1 5 -1;0 -1 0];
    Image = double(Image);
    [Height,Width] = size(Image);
    OutImage = zeros(Height,Width);
    for i = 2:Height-1
        for j = 2:Width-1
            OutImage(i,j) = abs(sum(sum(Image(i-1:i+1,j-1:j+1).*w1)));
        end
    end
end

边缘检测： 

clear;
clc;
%% laplace算子函数
% 主程序入口
Main();
%%
% 主函数
function Main()
    Image = imread('lena.tif');
    OutImage = uint8(Method_Laplace(Image));
    figure('name','Laplace边缘检测');
    imshow(OutImage);
    title('Laplace边缘');
end
%%
% laplace算子函数
function OutImage = Method_Laplace(Image)
    w1 = [0 1 0;1 -4 1;0 1 0];
    Image = double(Image);
    [Height,Width] = size(Image);
    OutImage = zeros(Height,Width);
    for i = 2:Height-1
        for j = 2:Width-1
            OutImage(i,j) = abs(sum(sum(Image(i-1:i+1,j-1:j+1).*w1)));
        end
    end
end

Robert算子与Sobel算子在此方面与之也是大同小异，仅仅是算子针对的差分方向不同而已，此处不再赘述。