空间域滤波：图像平滑和锐化

（1）自定义一个空间域平滑滤波函数，以达到滤除指定类型噪声(如高斯噪声和椒盐噪声等)的目的，该函数的输入参数包括滤波器类型filter_type (如高斯均值滤波、中值滤波、最大/小值滤波等)、处理前的图像ima(实验图像文件名为EXP4-1.jpg和EXP4-2.jpg)以及滤波器尺寸filter_size，输出参数为滤波处理后的图像im2。

EXP4-1.jpg

EXP4-2.jpg

自定义的滤波函数：

定义一个空间域平滑滤波函数 mySpatialFilter，该函数包含滤波器类型有3x3的高斯均值滤波器、中值滤波器、最大值滤波器以及最小值滤波器。

%pho 表示输入图像的名称
%filter_type 所选择的滤波器类型。有：高斯均值滤波器、中值滤波器、最大值滤波器、最小值滤波器
%filter_size 定义滤波器大小，如 filter_size=[3 3];
function ima2 = mySpatialFilter(pho, filter_type, filter_size)
ima=imread(pho); %根据路径读取原图像
ima=double(ima); %转化为双精度类型
[r c]=size(ima); %读取原图像的行数和列数
 
if strcmp(filter_type, '高斯均值滤波器')
    filter_size=[3 3];
end
m=filter_size(1); %读取设置的滤波器模板的大小
n=filter_size(2);
 
rs=r+m-1; %计算边界填充后的图像大小（填充模板大小）
cs=c+n-1;
tIma=zeros(rs,cs); %创建一个边界填充后的矩阵
tIma((m-1)/2+1:rs-(m-1)/2, (n-1)/2+1:cs-(n-1)/2)=ima; %矩阵中间填充原图像灰度级
 
%以最邻近的边界作为边界填充的灰度值
%先填充上边界
for i=1:(m-1)/2
    tIma(i,(n-1)/2+1:cs-(n-1)/2)=ima(1);
end
%填充下边界
for i=rs-(m-1)/2+1:rs
    tIma(i,(n-1)/2+1:cs-(n-1)/2)=ima(r);
end
%填充左边界
for i=1:(n-1)/2
    tIma(:,i)=tIma(:,(n-1)/2+1);
end
%填充右边界
for i=cs-(n-1)/2+1:cs
    tIma(:,i)=tIma(:,cs-(n-1)/2);
end
 
rIma=zeros(rs,cs);
%遍历该矩阵，套用滤波模板
for i=1:rs-m+1
    for j=1:cs-n+1
        %读取模板覆盖处的图像
        tempIma=tIma(i:i+m-1, j:j+n-1);
        %计算该模板覆盖的中心位置
        centerR=i+(m-1)/2;
        centerC=j+(m-1)/2;
        switch filter_type
            case '高斯均值滤波器'
                %定义3 x 3的高斯均值滤波器模板
                module=[1 2 1; 2 4 2; 1 2 1];
                module=module(:)';
                %为邻域以该模板为系数进行加权求平均，得出的值赋值得邻域中心
                rIma(centerR,centerC)=module*tempIma(:)/sum(module);
            case '中值滤波器'
                %为邻域内的值排序，求得中值，将该值赋值给邻域中心
                tempIma=sort(tempIma(:));
                rIma(centerR,centerC)=tempIma((length(tempIma)+1)/2,1);
            case '最大值滤波器'
                %求邻域中的最大灰度值，将该值赋值给邻域中心
                rIma(centerR,centerC)=max(tempIma(:));
            case '最小值滤波器'
                %求邻域中的最小灰度值，将该值赋值给邻域中心
                rIma(centerR,centerC)=min(tempIma(:));
            otherwise
                error('不存在该滤波器');
        end
    end
end
%去除原先填充的边界，得出最终结果
ima2=rIma((m-1)/2+1:rs-(m-1)/2, (n-1)/2+1:cs-(n-1)/2);
ima2=uint8(ima2);

定义一个使用各种滤波器对输入图像 EXP4-1.jpg 进行处理的函数 useMySpatialFilter1，输出结果。

function useMySpatialFilter1()
%读取 EXP4-1.jpg图像，进行显示
ima1=imread('EXP4-1.jpg');
subplot(151);
imshow(ima1);
title('输入图像 EXP4-1.jpg');
 
filter_size=[3 3]; %定义滤波器的大小为 3x3
type={'高斯均值滤波器' '中值滤波器' '最大值滤波器' '最小值滤波器'};
tit={'3x3 高斯均值滤波器' '3x3 中值滤波器' '3x3 最大值滤波器' '3x3 最小值滤波器'};
%运用各种滤波器，输出结果进行显示
for i=1:4
    ima2=mySpatialFilter('EXP4-1.jpg',type{i}, filter_size);
    subplot(1,5,i+1);
    imshow(ima2);
    title(tit(i));
end

输出结果如下：

结果分析：从输入图像中我们可以看出输入图像受椒盐噪声的干扰，那么在大小3x3的邻域中，椒盐噪声体现为邻域中的最大灰度值和最小灰度值。使用中值滤波器可以滤除最大灰度值和最小灰度值，以达到去除椒盐噪声的结果。这时使用高斯均值滤波器的效果就没有比中值滤波器的效果明显。因为高斯均值使用的是不同系数加权，即可能出现椒盐噪声处处于邻域中的最大权重，使得椒盐噪声有所加强。而最大值滤波器突出显示了邻域中灰度值最大部分，尤其突出了盐粒噪声，出现更多白点。最小值滤波器突出显示了邻域中灰度值最小部分，尤其突出了胡椒噪声，出现更多黑点。

定义一个使用各种滤波器对输入图像 EXP4-2.jpg 进行处理的函数 useMySpatialFilter2，输出结果。

function useMySpatialFilter2()
%读取 EXP4-2.jpg图像，进行显示
ima3=imread('EXP4-2.jpg');
subplot(151);
imshow(ima3);
title('输入图像 EXP4-2.jpg');
 
filter_size=[3 3];%定义滤波器的大小为 3x3
type={'高斯均值滤波器' '中值滤波器' '最大值滤波器' '最小值滤波器'};
tit={'3x3 高斯均值滤波器' '3x3 中值滤波器' '3x3 最大值滤波器' '3x3 最小值滤波器'};
%运用各种滤波器，输出结果进行显示
for i=1:4
    ima2=mySpatialFilter('EXP4-2.jpg',type{i}, filter_size);
    subplot(1,5,i+1);
    imshow(ima2);
    title(tit(i));
end

输出结果如下：

结果分析：输入图像受高斯噪声的影响，此时使用高斯均值滤波器和使用中值滤波器的处理效果相差不大。

（2）自定义一个空间域微分滤波器，要求该滤波器采用Sobel梯度算子获取输入图像的梯度图像，该滤波函数的输入参数为处理前的图像ima(实验图像文件名为EXP4-3.tif)，输出参数为输入图像的梯度图像im2。

定义一个空间域微分滤波器，使用Sobel梯度算子获取输入图像的梯度图像。

mySobelGradCal

%mySobelGradCal函数使用Sobel梯度算子作为微分滤波器模板，对输入图像进行锐化操作
%参数 pho 表示输入图像的路径，如 EXP4-3.tif
function grad2=mySobelGradCal(pho)
ima=imread(pho); %读取输入图像
ima=double(ima); %转化为双精度类型对图像进行操作
[r c]=size(ima); %获取输入的图像的行数和列数
m=3; %Sobel梯度算子的梯度模板大小（行）
n=3; %Sobel梯度算子的梯度模板大小（列）
%边界填充后的大小
rs=r+m-1;
cs=c+n-1;
%创建边界填充的矩阵
tIma=zeros(rs,cs);
%将输入图像填充在该矩阵的中心
tIma((m-1)/2+1:rs-(m-1)/2,(n-1)/2+1:cs-(n-1)/2)=ima;
 
for i=1:(m-1)/2
    %为矩阵填充下边界，根据输入图像的下边界进行填充
    tIma(i,(n-1)/2+1:cs-(n-1)/2)=ima(1,:);
    %为矩阵填充上边界，根据输入图像的上边界进行填充
    tIma(rs-i-1,(n-1)/2+1:cs-(n-1)/2)=ima(r,:);
end
 
for i=1:(n-1)/2
    %为矩阵填充左边界，根据矩阵已填充的左边一列进行填充
    tIma(:,i)=tIma(:,(n-1)/2+1);
    %为矩阵填充右边界，根据矩阵已填充的右边一列进行填充
    tIma(:,cs-i-1)=tIma(:,cs-(n-1)/2);
end
 
grad2=zeros(rs,cs);
%以可以让Sobel梯度模板覆盖tIma对tIma进行遍历
for i=1:rs-m+1
    for j=1:cs-n+1
        %取出当前被梯度模板覆盖的区域（灰度值）——邻域
        tempIma=tIma(i:i+m-1,j:j+n-1);
        %求邻域的中心
        centerR=i+(m-1)/2;
        centerC=j+(n-1)/2;
        %使用Sobel梯度算子进行计算
        x=abs(sum([tempIma(3,:) tempIma(3,2)])-sum([tempIma(1,:) tempIma(1,2)]));
        y=abs(sum([tempIma(:,3); tempIma(2,3)])-sum([tempIma(:,1); tempIma(2,1)]));
        grad2(centerR,centerC)=x+y;
    end
end
%截取出grad2中中心部分即为输出结果
grad2=grad2((m-1)/2+1:rs-(m-1)/2,(n-1)/2+1:cs-(n-1)/2);
grad2=uint8(grad2);

useMySobelGradCal

%useMySobelGradCal函数显示输入图像和对输入图像应用mySobelGradCal函数后所得图像
%参数pho 表示输入图像的路径，如 EXP4-3.tif
function []=useMySobelGradCal(pho)
%读取输入图像
ima1=imread(pho);
%显示输入图像
subplot(121);
imshow(ima1);
title('输入图像');
 
%对输入图像应用mySobelGradCal函数得到Sobel梯度图像
ima2=mySobelGradCal(pho);
%显示Sobel梯度图像
subplot(122);
imshow(ima2);
title('Sobel梯度图像');

输出结果如下：

需要注意的是：上述的输入图像是受噪声污染的图像，所以在使用微分滤波器前，可以先使用此前设计的滤波器合理地处理一下图像，再进行微分滤波器的处理。