求一个用均值滤波,中值滤波和维纳滤波分别去除高斯噪声和椒盐噪声的 matlab程序
psnr:
function [PSNR, MSE] = psnr(X, Y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 计算峰值信噪比PSNR
% 将RGB转成YCbCr格式进行计算
% 如果直接计算会比转后计算值要小2dB左右(当然是个别测试)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(X,3)~=1 %判断图像时不是彩色图,如果是,结果为3,否则为1
org=rgb2ycbcr(X);
test=rgb2ycbcr(Y);
Y1=org(:,:,1);
Y2=test(:,:,1);
Y1=double(Y1); %计算平方时候需要转成double类型,否则uchar类型会丢失数据
Y2=double(Y2);
else %灰度图像,不用转换
Y1=double(X);
Y2=double(Y);
end
if nargin<2
D = Y1;
else
if any(size(Y1)~=size(Y2))
error('The input size is not equal to each other!');
end
D = Y1 - Y2;
end
MSE = sum(D(:).*D(:)) / numel(Y1);
PSNR = 10*log10(255^2 / MSE);
你可以参考下此文章matlab实现PSNR_御坂御坂001的博客-CSDN博客
a=imread('mengnalisha.jpg');
b=rgb2gray(a);
subplot(241),imshow(a),title('原图');
subplot(242),imshow(b),title('灰度图');
c=imnoise(b,'salt & pepper',0.02);
subplot(243),imshow(c),title('添加椒盐噪声图像');
x=ones(3,3)/9; %定义一个3*3的均值滤波器
d=imfilter(c,x); %滤波
subplot(244),imshow(d),title('均值滤波后的图像');
e=medfilt2(c); %中值滤波
subplot(245),imshow(e),title('中值滤波后的图像');
f=wiener2(c,[3 3]); %维纳滤波,[3 3]是设置的邻域大小
subplot(246),imshow(f),title('维纳滤波后的图像');
[C1,S1] = wavedec2(c,2,'coif3'); %提取小波分解中第一层的低频图像,即实现了低通滤波去噪
n = [1,2]; %设置尺度向量n
p = [10.12,23.28]; %设置阈值向量p
nc1 = wthcoef2('h',C1,S1,n,p,'s');%对三个方向h,v,d(水平、垂直、对角线)高频系数进行软阈值处理
nc1 = wthcoef2('v',nc1,S1,n,p,'s');
nc1 = wthcoef2('d',nc1,S1,n,p,'s');
x1 = waverec2(nc1,S1,'coif3'); %对新的小波分解结构[c,s]进行重构
Z1=uint8(x1); %将double型改为unit8型显示图像
subplot(247),imshow(Z1),title('小波软阈值处理后的图像');
nc2 = wthcoef2('h',C1,S1,n,p,'h');%三个方向硬阈值处理
nc2 = wthcoef2('v',nc2,S1,n,p,'h');
nc2= wthcoef2('d',nc2,S1,n,p,'h');
x2 = waverec2(nc2,S1,'coif3');
Z2=uint8(x2);
subplot(248),imshow(Z2),title('小波硬阈值处理后的图像');
[a1,a2]=psnr(Z1,b); %求软阈值去噪的峰值信噪比a1和信噪比a2
[a3,a4]=psnr(Z2,b); %求硬阈值去噪的峰值信噪比a1和信噪比a2
% RMSE=sqrt((sum((c-b).^2))./2);
%软阈值结果评价
B=double(b);
A=B-x1;
MSE = sum(A(:).*A(:))/numel(B); %均方误差MSE,numel()函数返回矩阵元素个数
RMSE=sqrt(MSE);
SNR = 10*log10(sum(B(:).*B(:))/MSE/numel(B));%信噪比SNR
PSNR = 10*log10(255^2/MSE); %峰值信噪比PSNR
%硬阈值结果评价
C=B-x2;
MSE2 = sum(C(:).*C(:))/numel(B);
RMSE2=sqrt(MSE2);
SNR2 = 10*log10(sum(B(:).*B(:))/MSE2/numel(B));
PSNR2 = 10*log10(255^2/MSE2);
%均值滤波结果评价
D=double(d);
D1=B-D;
MSE3 = sum(D1(:).*D1(:))/numel(B);
RMSE3=sqrt(MSE3);
SNR3 = 10*log10(sum(B(:).*B(:))/MSE3/numel(B));
PSNR3 = 10*log10(255^2/MSE3);
%中值滤波结果评价
E=double(e);
E1=B-E;
MSE4 = sum(E1(:).*E1(:))/numel(B);
RMSE4=sqrt(MSE4);
SNR4 = 10*log10(sum(B(:).*B(:))/MSE4/numel(B));
PSNR4 = 10*log10(255^2/MSE4);
%维纳滤波结果评价
F=double(f);
F1=B-F;
MSE5 = sum(F1(:).*F1(:))/numel(B);
RMSE5=sqrt(MSE5);
SNR5 = 10*log10(sum(B(:).*B(:))/MSE5/numel(B));
PSNR5 = 10*log10(255^2/MSE5);
如果去除高斯噪声的将原图加入高斯噪声:
c=imnoise(b,'gaussian',0.02);
subplot(243),imshow(c),title('添加高斯噪声图像');
再去噪
您好,我是有问必答小助手,您的问题已经有小伙伴解答了,您看下是否解决,可以追评进行沟通哦~
如果有您比较满意的答案 / 帮您提供解决思路的答案,可以点击【采纳】按钮,给回答的小伙伴一些鼓励哦~~
ps:问答VIP仅需29元,即可享受5次/月 有问必答服务,了解详情>>>https://vip.csdn.net/askvip?utm_source=1146287632