[數字影象處理]影象去噪初步(1)--均值濾波器
阿新 • • 發佈：2019-02-10
close all; clear all; clc; f = imread('./ckt-board-orig.tif'); f = mat2gray(f,[0 255]); [M,N] = size(f); %% ---------------gaussian noise------------------- a = 0; b = 0.1; n_gaussian = a + b .* randn(M,N); g_gaussian = f + n_gaussian; g_gaussian(find(g_gaussian > 1)) = 1; g_gaussian(find(g_gaussian < 0)) = 0; figure(); subplot(1,2,1); imshow(g_gaussian,[0 1]); xlabel('a).Image corrupted by Gaussian noise'); subplot(1,2,2); x = linspace(-0.2,1.2,358); h = hist(g_gaussian,x)/(M*N); Histogram = zeros(358,1); for y = 1:256 Histogram = Histogram + h(:,y); end bar(-0.2:1/255:1.2,Histogram); axis([-0.2 1.2 0 0.014]),grid; xlabel('b).The Histogram of a'); ylabel('Number of pixels'); g_diff = abs(g_gaussian - f); MSE = sum(sum(g_diff .^2))/(M*N); PSNR = 10*log10((1*1)/MSE) [SSIM_G mssim] = ssim_index(f,g_gaussian,[0.01 0.03],ones(8),1); %% ---------------Denoise(Gaussian noise)------------------- g_Ex = zeros(M+2,N+2); for x = 1:M g_Ex(x+1,:) = [g_gaussian(x,1) g_gaussian(x,:) g_gaussian(x,N)]; end g_Ex(1,:) = g_Ex(2,:); g_Ex(M+2,:) = g_Ex(M+1,:); % Arithemtic mean filter g_amf = zeros(M,N); for x = 2:M+1 for y = 2:N+1 g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex(x ,y); g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex(x-1,y); g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex(x ,y-1); g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex(x+1,y); g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex( x,y+1); g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex(x-1,y+1); g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex(x+1,y-1); g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex(x+1,y+1); g_amf(x-1,y-1) = g_amf(x-1,y-1) + g_Ex(x-1,y-1); g_amf(x-1,y-1) = g_amf(x-1,y-1)/9; end end g_amf_diff = abs(g_amf - f); MSE_amf = sum(sum(g_amf_diff .^2))/(M*N); PSNR_amf = 10*log10((1*1)/MSE_amf) [SSIM_amf mssim] = ssim_index(f,g_amf ,[0.01 0.03],ones(8),1); % Geometric mean filter g_gmf = zeros(M,N); for x = 2:M+1 for y = 2:N+1 g_gmf(x-1,y-1) = g_Ex(x ,y); g_gmf(x-1,y-1) = g_gmf(x-1,y-1) * g_Ex(x-1,y); g_gmf(x-1,y-1) = g_gmf(x-1,y-1) * g_Ex(x ,y-1); g_gmf(x-1,y-1) = g_gmf(x-1,y-1) * g_Ex(x+1,y); g_gmf(x-1,y-1) = g_gmf(x-1,y-1) * g_Ex( x,y+1); g_gmf(x-1,y-1) = g_gmf(x-1,y-1) * g_Ex(x-1,y+1); g_gmf(x-1,y-1) = g_gmf(x-1,y-1) * g_Ex(x+1,y-1); g_gmf(x-1,y-1) = g_gmf(x-1,y-1) * g_Ex(x+1,y+1); g_gmf(x-1,y-1) = g_gmf(x-1,y-1) * g_Ex(x-1,y-1); end end g_gmf = (g_gmf).^(1/9); g_gmf_diff = abs(g_gmf - f); MSE_gmf = sum(sum(g_gmf_diff .^2))/(M*N); PSNR_gmf = 10*log10((1*1)/MSE_gmf) [SSIM_gmf mssim] = ssim_index(f,g_gmf ,[0.01 0.03],ones(8),1); figure(); subplot(1,2,1); imshow(g_amf,[0 1]); xlabel('a).Ruselt of Denoise by Amf'); figure(); subplot(1,2,2); imshow(g_gmf,[0 1]); xlabel('a).Ruselt of Denoise by Gmf'); % Harmonic mean filter g_hmf = zeros(M,N); for x = 2:M+1 for y = 2:N+1 g_hmf(x-1,y-1) = 1/g_Ex(x ,y); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x-1,y); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x ,y-1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x+1,y); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex( x,y+1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x-1,y+1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x+1,y-1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x+1,y+1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x-1,y-1); g_hmf(x-1,y-1) = 9/g_hmf(x-1,y-1); end end g_hmf_diff = abs(g_hmf - f); MSE_hmf = sum(sum(g_hmf_diff .^2))/(M*N); PSNR_hmf = 10*log10((1*1)/MSE_hmf) [SSIM_hmf mssim] = ssim_index(f,g_hmf ,[0.01 0.03],ones(8),1); figure(); subplot(1,2,1); imshow(g_hmf,[0 1]); xlabel('c).Ruselt of Denoise by Hmf'); %% ---------------Denoise(salt and pepper noise)---------- close all; clear all; clc; f = imread('./ckt-board-orig.tif'); f = mat2gray(f,[0 255]); [M,N] = size(f); %% ---------------Salt & pepper------------------- b = 0; %salt a = 0.2; %pepper x = rand(M,N); g_sp = zeros(M,N); g_sp = f; g_sp(find(x<=a)) = 0; g_sp(find(x > a & x<(a+b))) = 1; g_diff = abs(g_sp - f); MSE = sum(sum(g_diff .^2))/(M*N); PSNR = 10*log10((1*1)/MSE) figure(); subplot(1,2,1); imshow(g_sp,[0 1]); xlabel('a).Image corrupted by salt&pepper noise'); %% ---------------Denoise(Salt & pepper)------------------- g_Ex = zeros(M+2,N+2); for x = 1:M g_Ex(x+1,:) = [g_sp(x,1) g_sp(x,:) g_sp(x,N)]; end g_Ex(1,:) = g_Ex(2,:); g_Ex(M+2,:) = g_Ex(M+1,:); % Harmonic mean filter g_hmf = zeros(M,N); for x = 2:M+1 for y = 2:N+1 g_hmf(x-1,y-1) = 1/g_Ex(x ,y); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x-1,y); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x ,y-1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x+1,y); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex( x,y+1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x-1,y+1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x+1,y-1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x+1,y+1); g_hmf(x-1,y-1) = g_hmf(x-1,y-1) + 1/g_Ex(x-1,y-1); g_hmf(x-1,y-1) = 9/g_hmf(x-1,y-1); end end g_hmf_diff = abs(g_hmf - f); MSE_hmf = sum(sum(g_hmf_diff .^2))/(M*N); PSNR_hmf = 10*log10((1*1)/MSE_hmf) figure(); subplot(1,2,1); imshow(g_sp,[0 1]); xlabel('a).Image corrupted by pepper noise'); subplot(1,2,2); imshow(g_hmf,[0 1]); xlabel('b).Ruselt of Denoise by Hmf'); % Contraharmonic mean filter Q = -1.5; g_cmf = zeros(M,N); for x = 2:M+1 for y = 2:N+1 g_cmf_M = (g_Ex(x ,y))^(1+Q); g_cmf_M = g_cmf_M + (g_Ex(x-1,y))^(1+Q); g_cmf_M = g_cmf_M + (g_Ex(x ,y-1))^(1+Q); g_cmf_M = g_cmf_M + (g_Ex(x+1,y))^(1+Q); g_cmf_M = g_cmf_M + (g_Ex( x,y+1))^(1+Q); g_cmf_M = g_cmf_M + (g_Ex(x-1,y+1))^(1+Q); g_cmf_M = g_cmf_M + (g_Ex(x+1,y-1))^(1+Q); g_cmf_M = g_cmf_M + (g_Ex(x+1,y+1))^(1+Q); g_cmf_M = g_cmf_M + (g_Ex(x-1,y-1))^(1+Q); g_cmf_D = (g_Ex(x ,y))^(Q); g_cmf_D = g_cmf_D + (g_Ex(x-1,y))^(Q); g_cmf_D = g_cmf_D + (g_Ex(x ,y-1))^(Q); g_cmf_D = g_cmf_D + (g_Ex(x+1,y))^(Q); g_cmf_D = g_cmf_D + (g_Ex( x,y+1))^(Q); g_cmf_D = g_cmf_D + (g_Ex(x-1,y+1))^(Q); g_cmf_D = g_cmf_D + (g_Ex(x+1,y-1))^(Q); g_cmf_D = g_cmf_D + (g_Ex(x+1,y+1))^(Q); g_cmf_D = g_cmf_D + (g_Ex(x-1,y-1))^(Q); g_cmf(x-1,y-1) = g_cmf_M/g_cmf_D; end end g_cmf_diff = abs(g_cmf - f); MSE_cmf = sum(sum(g_cmf_diff .^2))/(M*N); PSNR_cmf = 10*log10((1*1)/MSE_cmf) figure(); subplot(1,2,1); imshow(g_cmf,[0 1]); xlabel('b).Ruselt of Denoise by Cmf(Q=-1.5)');