clear
clc

X1=imread('Pic1_3.bmp');

figure
subplot(131)
imshow(X1)
subplot(132)
imhist(X1)
XX1=im2bw(X1,254/255);
subplot(133)
for kk=1:2
m=30;
for i=m+1:size(XX1,1)-m
    for j=m+1:size(XX1,2)-m
        a=XX1(i-m:i+m,j-m:j+m);
        if length(find(a(1:m,m+1)==0))>0 & length(find(a(m+1:end,m+1)==0))>0 & length(find(a(m+1,1:m)==0))>0 & length(find(a(m+1,m+1:end)==0))>0
            XX1(i,j)=0;
        end
    end
end
end
imshow(XX1)

Ig=XX1;
%获取算子
s=GetStrelList();
%串联去噪
e=ErodeList(Ig,s);
%计算权重
f=GetRateList(Ig,e);
%并联
Igo=GetRemoveResult(f,e);

%计算PSNR值
psnr1=PSNR(XX1,e.eroded_co12);
psnr2=PSNR(XX1,e.eroded_co22);
psnr3=PSNR(XX1,e.eroded_co32);
psnr4=PSNR(XX1,e.eroded_co42);
psnr5=PSNR(XX1,Igo);
psnr_list=[psnr1 psnr2 psnr3 psnr4 psnr5];
M{1,1}=e.eroded_co12;
M{1,2}=e.eroded_co22;
M{1,3}=e.eroded_co32;
M{1,4}=e.eroded_co42;
M{1,5}=Igo;

figure
plot(1:5,psnr_list,'r+-');
set(gca,'XTick',0:6,'XTickLabel',{'','串联1','串联2','串联3','串联4','并联',''});
grid on;
title('PSNR曲线比较')

[~,b]=max(psnr_list);
result=M{1,b};


figure

imshow(result)

imhist(result)
result=im2bw(result,0.3);

Y1=bwperim(result);
imshow(Y1)
imwrite(Y1,'pic1.bmp')

function s=GetStrelList()
%获取算子
%输出参数
%s—算子结构体

%生成串联算子
s.co11=strel('line',3,-45);
s.co12=strel('line',5,-45);
%生成串联算子
s.co21=strel('line',3,45);
s.co22=strel('line',5,45);
%生成串联算子
s.co31=strel('line',3,90);
s.co32=strel('line',5,90);
%生成串联算子
s.co41=strel('line',3,0);
s.co42=strel('line',5,0);

function Igo=GetRemoveResult(f,e)
%并联去噪
%输入参数
% f—权值向量
% e—串联结果
%输出参数
% Igo—处理结果
Igo=...
    f.df1/f.df*double(e.eroded_co12)+f.df2/f.df*double(e.eroded_co22)+...
    f.df3/f.df*double(e.eroded_co32)+f.df4/f.df*double(e.eroded_co42);
Igo=mat2gray(Igo);

function f=GetRateList(Ig,e)
%计算权重
%输入参数
% Ig—图像矩阵
% e—串联结果
%输出参数
% f—处理结果
f.df1=sum(sum(abs(double(e.eroded_co12)-double(Ig))));
f.df2=sum(sum(abs(double(e.eroded_co22)-double(Ig))));
f.df3=sum(sum(abs(double(e.eroded_co32)-double(Ig))));
f.df4=sum(sum(abs(double(e.eroded_co42)-double(Ig))));
f.df=sum([f.df1 f.df2 f.df3 f.df4]);

function e=ErodeList(Ig,s)
%串联去噪
%输入参数
% Ig—图像矩阵
% s—算子
%输出参数
% e—处理结果
e.eroded_co11=imerode(Ig,s.co11);
e.eroded_co12=imerode(e.eroded_co11,s.co12);
e.eroded_co21=imerode(Ig,s.co21);
e.eroded_co22=imerode(e.eroded_co21,s.co22);
e.eroded_co31=imerode(Ig,s.co31);
e.eroded_co32=imerode(e.eroded_co21,s.co32);
e.eroded_co41=imerode(Ig,s.co41);
e.eroded_co42=imerode(e.eroded_co21,s.co42);

function S=PSNR(s,t)
%计算PSNR
%输入参数
% S—图像矩阵1
% t—图像矩阵2
%输出参数
% S—结果

%预处理
[m,n,~]=size(s);
s=im2uint8(mat2gray(s));
t=im2uint8(mat2gray(t));
s=double(s);
t=double(t);
%初值
sd=0;
mi=m*n*max(max(s.^2));
%计算
for u=1:m
    for v=1:n
        sd=sd+(s(u,v)-t(u,v))^2;
    end
end
if sd==0
    sd=1;
end
S=mi/sd;
S=10*log10(S);