【数字图像处理】Spatial Filtering and Median Filtering

实验要求：

编写一个能够完成两幅图像之间加、减、乘、除四种算术运算的函数。另外，对于两幅图像的乘法，所编写的乘法程序还要能够完成一幅图像乘以一个常数的功能。使用图 Fig1.10(4)和 Fig1.10(5)验证实验。

编写一个对图像进行空域滤波的函数(实现方法见课本 3.5 节)。空域模板的尺寸固定为 3 x 3 大小, 但模板
中的系数要做为程序的输入参数。

(1.a) 使用上述函数实现图像的平滑滤波，其中的模板使用图 3.32 所示的掩模。
(1.b) 使用上述函数实现公式(3.6-8)和公式(3.6-9)描述的拉普拉斯图像锐化增强技术，其中的拉普拉斯模板
使用图 3.37(c)(d)所示的掩模。
(1.c) 使用上述函数实现公式(3.6-8)和公式(3.6-9)描述的一阶导数锐化增强技术，其中的一阶导数模板使用
图 3.41(d)(e)所示的掩模。
(1.d) 使用课本图 3.38(a)，对(1.a)(1.b)(1.c)中的程序实现进行验证,并对结果进行简要对比分析。

(2.a) 编写一个给图像中添加椒盐噪声的程序，程序的输入参数为两个噪声分量的概率值。
(2.b) 编写的程序，实现 3 x 3 中值滤波。
(2.c) 给图 5.7(a)添加 Pa = Pb = 0.2 的椒盐噪声。
(2.d) 对 (2.c)中得到的图像进行中值滤波。说明实验结果与课本图 5.10(b)的主要差异。

实现代码：

1.图像间四则算术运算部分的程序：
% 主函数Image_arithmetic_test.m
% 说明：编写一个能够完成两幅图像之间加、减、乘、除四种算术运算的函数。另外，对于两幅图像的乘法，所编写的乘法程序还要能够完成一幅图像乘以一个常数的功能。使用图 Fig1.10(4)和 Fig1.10(5)验证实验。
clear all; close all; clc;
 
img_A = imread('Fig1.10(4).jpg');
img_B = imread('Fig1.10(5).jpg');
 
img1 = compute(img_A,'+',img_B);                     % 加法
img2 = compute(img_A,'-',img_B);                     % 减法
img3 = compute(img_A,'*',img_B);                     % 乘法
img4 = compute(img_A,'/',img_B);                     % 除法
img5 = compute(img_A,'*',0.5,1);                     % 图像乘以0.5
img6 = compute(img_A,'*',2,1);                       % 图像乘以2
 
figure;
subplot(2,4,1);imshow(img5);title('(a)图像Fig1.10(4)');
subplot(2,4,2);imshow(img6);title('(b)图像Fig1.10(5)');
subplot(2,4,3);imshow(img1);title('(c)图像相加');
subplot(2,4,4);imshow(img2);title('(d)图像相减');
subplot(2,4,5);imshow(img3);title('(e)图像相乘');
subplot(2,4,6);imshow(img4);title('(f)图像相除');
subplot(2,4,7);imshow(img5);title('(g)图像Fig1.10(4)乘以0.5');
subplot(2,4,8);imshow(img6);title('(h)图像Fig1.10(4)乘以2');

function img_out = compute(A,param,B,is_num)
 
% 该函数用于完成两幅图像之间的加、减、乘、除四种算术运算，
% 其中乘法还能实现一幅图像与一个常数相乘。
%  A为输入的第一幅图像，B为输入的第二幅图像或者常数，由is_num确定，
%  param指+、-、*、/四个运算符。
[M,N]= size(A);
img_out = zeros(M,N);
 
if nargin == 3 
    for m = 1:M
        for n = 1:N   
            if param == '+'                          % 加法
                img_out(m,n) = A(m,n)+B(m,n);
            elseif param == '-'                      % 减法
                img_out(m,n) = abs(A(m,n)-B(m,n));
            elseif param == '*'                      % 乘法
                img_out(m,n) = A(m,n)*B(m,n);                
            elseif param == '/'                      % 除法(考虑除数是否为零)
                if(B(m,n) ~= 0)                         
                    img_out(m,n) = A(m,n)/B(m,n);
                else
                    img_out(m,n) = 255;
                end
            end
        end
    end    
elseif nargin ==4 && param == '*' && is_num == 1     % 图像乘以一个常数
    for m = 1:M
        for n = 1:N
            img_out(m,n) = A(m,n)*B; 
        end
    end    
else
    disp('Input Error');
end
img_out = uint8(img_out);
end

2.图像空域滤波部分的程序：
% 主函数Smooth_filter_test.m
%说明：编写一个对图像进行空域滤波的函数(实现方法见课本 3.5 节)。空域模板的尺寸固定为 3 x 3 大小, 但模板中的系数要做为程序的输入参数。
clear all; close all; clc;
 
img = imread('Fig3.38(a).jpg');
 
M_smooth1 = [1 1 1;1 1 1;1 1 1]/9;                    % 平滑滤波的掩模
M_smooth2 = [1 2 1;2 4 2;1 2 1]/16;
M_lap1 = [0 -1 0;-1 -4 -1;0 -1 0];                    % 拉普拉斯推向锐化增强技术的掩模
M_lap2 = [-1 -1 -1;-1 8 -1;-1 -1 -1];
M_sobel1 = [-1 -2 -1;0 0 0;1 2 1];                    % 一阶导数锐化增强技术的掩模
M_sobel2 = [-1 0 1;-2 0 2;-1 0 1];               
 
img_smooth1 = filter_spatial(img,M_smooth1);          % 平滑滤波
img_smooth2 = filter_spatial(img,M_smooth2);
 
f_lap1 = filter_spatial(img,M_lap1);                  % 拉普拉斯图像锐化增强
mask_lap1 = compute(img,'-',f_lap1);
img_lap1 = compute(img,'+', mask_lap1);
 
f_lap2 = filter_spatial(img,M_lap2);
mask_lap2 = compute(img,'-',f_lap2);
img_lap2 = compute(img,'+', mask_lap2);
 
f_sobel1 = filter_spatial(img,M_sobel1);              % 一阶导数锐化增强
mask_sobel1 = compute(img,'-',f_sobel1);
img_sobel1 = compute(img,'+', mask_sobel1);
 
f_sobel2 = filter_spatial(img,M_sobel2);
mask_sobel2 = compute(img,'-',f_sobel2);
img_sobel2 = compute(img,'+', mask_sobel2);
 
subplot(3,3,1);imshow(img);title('(a)：原图像Fig3.38(a)');      
subplot(3,3,4);imshow(img_smooth1);title('(b)：平滑滤波：模板a');
subplot(3,3,5);imshow(img_smooth2);title('(c)：平滑滤波：模板b');
subplot(3,3,6);imshow(img_lap1);title('(d)：拉普拉斯锐化：模板c');
subplot(3,3,7);imshow(img_lap2);title('(e)：拉普拉斯锐化：模板d');
subplot(3,3,8);imshow(img_sobel1);title('(f)：sobel锐化：模板d');
subplot(3,3,9);imshow(img_sobel2);title('(g)：sobel锐化：模板e');

% 该函数用于拉普拉斯图像锐化增强，给Smooth_filter_test.m程序调用
function img_out = filter_spatial(img_in,M)
% 该函数用于对图像进行空间滤波，空域模板M的尺寸固定为3 x 3 大小。
[X,Y] = size(img_in);
img_out = zeros(X,Y);
 
for x = 2:X-1
    for y = 2:Y-1
        img_out(x,y) = img_in(x-1,y-1)*M(1,1)+img_in(x-1,y)*M(1,2)+img_in(x-1,y+1)*M(1,3)...
            +img_in(x,y-1)*M(2,1)+img_in(x,y)*M(2,2)+img_in(x,y+1)*M(2,3)...
            +img_in(x+1,y-1)*M(3,1)+img_in(x+1,y)*M(3,2)+img_in(x+1,y+1)*M(3,3);
    end
end
 
img_out = uint8(img_out);
end    

3.图像加椒盐噪声和中值滤波部分的程序：
% 主函数Median_filtering_test.m
clear all; close all; clc;
 
img = imread('Fig5.07(a).jpg');
img_sp = salt_pepper(img,0.2,0.2);          % 椒盐噪声
img_med = filter_median(img_sp);            % 中值滤波
 
subplot(1,3,1);imshow(img);title('原图像Fig5.07(a)');
subplot(1,3,2);imshow(img_sp);title('添加椒盐噪声');
subplot(1,3,3);imshow(img_med);title('中值滤波');

function img_out = salt_pepper(img_in,Pa,Pb)
% 该函数用于给图像添加椒盐噪声。
 
if Pa>=1 || Pb>=1                   % 检测输入参数
    error('Input Error');
end
 
[M,N] = size(img_in);
img_out = img_in;
s = rand(M,N);                      % 产生M行N列的位于(0,1)区间的随机数
p = rand(M,N);
 
for m = 1:M
    for n = 1:N
        if s(m,n)<=Pa
            img_out(m,n)=0;         % 盐粒噪声
        end
        if p(m,n)<=Pb
            img_out(m,n)=255;       % 胡椒噪声
        end
    end
end
end

function img_out = filter_median(img_in)
% 该函数用于实现中值滤波
 
[R,C] = size(img_in);
img_in = double(img_in);
img_out = zeros(R,C);
 
for m = 2:R-1                   
    for n = 2:C-1
        M = [img_in(m-1,n-1),img_in(m-1,n),img_in(m-1,n+1);...
            img_in(m,n-1),img_in(m,n),img_in(m,n+1);...
            img_in(m+1,n-1),img_in(m+1,n),img_in(m+1,n+1)];
        img_out(m,n) = median(M);   %调用median函数进行冒泡排序
    end
end
 
img_out = uint8(img_out);
end        

function m = median(M)
% 该函数实现冒泡排序
% 该函数用于计算3×3矩阵的中值。
for i = 1:8
    for j = 1:8
        if M(j)>M(j+1)
            t = M(j);
            M(j) = M(j+1);
            M(j+1) = t;
        end
    end
end
 
m = M(5);
end