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

（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];endm=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    endend%去除原先填充的边界，得出最终结果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]; %定义滤波器的大小为 3x3type={'高斯均值滤波器' '中值滤波器' '最大值滤波器' '最小值滤波器'};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];%定义滤波器的大小为 3x3type={'高斯均值滤波器' '中值滤波器' '最大值滤波器' '最小值滤波器'};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.tiffunction 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;    endend%截取出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.tiffunction []=useMySobelGradCal(pho)%读取输入图像ima1=imread(pho);%显示输入图像subplot(121);imshow(ima1);title('输入图像'); %对输入图像应用mySobelGradCal函数得到Sobel梯度图像ima2=mySobelGradCal(pho);%显示Sobel梯度图像subplot(122);imshow(ima2);title('Sobel梯度图像');

输出结果如下：

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