《数字图像处理》冈萨雷斯学习笔记4：图像复原

g = imnoise(f,type,parameters)

g = imnoise(f,'gaussian',m,var) %将均值为m、方差为var的高斯噪声加到图像f上，默认值为均值是0，方差是0.01的噪声
g = imnoise(f,'salt & pepper',d) %用椒盐噪声污染图像f，其中d是噪声密度
g = imnoise(f,'poisson') %从数据中生成泊松噪声，而不是人工添加

g = imnoise2(type,M,N,a,b)

clc
clear
r1=imnoise2('gaussian',100000,1,1);%imnoise2是产生噪声矩阵，这里是产生高斯噪声，矩阵大小为10000*1
bins=100;                           %均值为0，方差为1
figure,hist(r1,bins);%将r矩阵中的数值直方图表示出来，共100个bin
title('guassian');
r2=imnoise2('uniform',1);%同上，此处产生的是0~1的均匀分布噪声，矩阵大小为100000*1
figure,hist(r2,bins);
title('uniform');
r3=imnoise2('salt & pepper',1000,0.1,0.27);%注意这里的salt & pepper中间记得留空格。                      
figure,hist(r3,bins);
title('salt & pepper');
r4=imnoise2('lognormal',1);%产生对数正态噪声，大小100000*1，偏移值默认为1，形状参数默认0.25
figure,hist(r4,bins);
title('lognormal');
r5=imnoise2('rayleigh',1);%rayleigh噪声，该噪声好像是对整体矢量合成的一种分布
figure,hist(r5,bins);
title('rayleigh');
r6=imnoise2('exponential',1);%指数分布的参数默认为1
figure,hist(r6,bins);
title('exponential');
% figure,subplot(321),imshow(r1),title('高斯噪声');
% subplot(322),imshow(r2),title('均匀分布噪声');
% subplot(323),imshow(r3),title('椒盐噪声');
% subplot(324),imshow(r4),title('对数正态噪声');
% subplot(325),imshow(r5),title('瑞利噪声');
% subplot(326),imshow(r6),title('指数噪声');

%% 空间噪声滤波
clc
clear
f=imread('.imagesdipum_images_ch05Fig0507(a)(checkeboard8_pixeldup_8).tif');
imshow(f),title('原始图像');

[M N]=size(f);
R=imnoise2('salt & pepper',0);%仅仅产生0.1的椒噪声
c=find(R==0);%把R中产生的椒点找出来
gp=f;
gp(c)=0;%在原始图像对应椒点的位置赋值为0
% figure,imshow(gp);
% title('被概率为0.1的胡椒噪声污染的图像');

R=imnoise2('salt & pepper',0.1);%仅仅产生0.1的盐噪声
c=find(R==1);%把R中产生的盐点找出来
gs=f;
gs(c)=255;%在原始图像对应盐点的位置赋值为1
% figure,imshow(gs);
% title('被概率为0.1的盐噪声污染的图像');

fp=spfilt(gp,'chmean',3,1.5);%反调和滤波，尺寸大小为3*3，Q默认为1.5,正的Q过滤胡椒噪声
% figure,imshow(fp);
% title('正Q反调和滤波器滤波的结果')


fs=spfilt(gs,-1.5);%负的Q过滤盐粒噪声
% figure,imshow(fs);
% title('负Q反调和滤波器滤波的结果')


fpmax=spfilt(gp,'max',3);
% figure,imshow(fpmax);
% title('最大值滤波后的结果');


fsmin=spfilt(gs,'min',imshow(fsmin);
% title('最小值滤波后的结果');
figure,imshow(gp);title('被概率为0.1的胡椒噪声污染的图像');
subplot(322),imshow(gs);title('被概率为0.1的盐噪声污染的图像');
subplot(323),imshow(fp);title('正Q反调和滤波器滤波的结果')
subplot(324),imshow(fs);title('负Q反调和滤波器滤波的结果')
subplot(325),imshow(fpmax);title('最大值滤波后的结果');
subplot(326),imshow(fsmin);title('最小值滤波后的结果');

f = adpmedian (g,Smax)

%% 自适应中值滤波器
clc
clear
f=imread('.imagesdipum_images_ch05Fig0507(a)(checkeboard8_pixeldup_8).tif');
figure,subplot(221),imshow(f),title('原始图像');

g=imnoise(f,0.25);%噪声点有白有黑，因为这是用的imnoise,不是imnoise2和imnoise3
subplot(222),imshow(g);                
title('被概率为0.25的椒盐噪声污染过的图像');

f1=medfilt2(g,[7 7],'symmetric'); %边缘扩展模式为对称扩展
subplot(223),imshow(f1);
title('用普通中值滤波器滤波后的图像');

f2=adpmedian(g,7);
subplot(224),imshow(f2);
title('用Smax=7的自适应中值滤波器滤波后的图像');

%% 模糊噪声图像建模
clc
clear
f=checkerboard(8);%直接生产黑白棋盘，每个小正方形一边的像素点为8个
figure,imshow(f);
title('原始图像');

PSF=fspecial('motion',7,45);%创建一个motion滤波器，长度为7，仰角为45度。
gb=imfilter(f,PSF,'circular');%具体motion滤波器是什么还不是很清楚。
subplot(222),imshow(gb);
title('使用motion滤波器模糊后');

noise=imnoise(zeros(size(f)),0.001);%产生高斯噪声，且后面要用到该噪声图像
subplot(223),imshow(noise,[]);
% subplot(223),imshow(noise); %全是黑的
title('高斯纯噪声')

g=gb+noise;
subplot(224),imshow(g,[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');

%==========================================================================
%%使用deconvwnr函数复原模糊噪声图像

figure,[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');
fr1=deconvwnr(g,PSF);%维纳滤波，去模糊化
subplot(222),imshow(fr1,[]);
title('简单维纳滤波后的结果');

Sn=abs(fft2(noise)).^2;%噪声功率谱
nA=sum(Sn(:))/prod(size(noise));%平均噪声功率谱，prod是元素相乘的意思，这里就是noise的长和宽相乘
Sf=abs(fft2(f)).^2;%信号功率谱
fA=sum(Sf(:))/prod(size(f));%平均信号功率谱
R=nA/fA;%求得信噪比

fr2=deconvwnr(g,R);%信噪比为常数的参数维纳滤波器
subplot(223),imshow(fr2,[]);
title('常数比率信噪比维纳滤波器');

NCORR=fftshift(real(ifft(Sn)));%噪声的自相关函数，why？
ICORR=fftshift(real(ifft(Sf)));%信号的自相关函数，why？
fr3=deconvwnr(g,NCORR,ICORR);%自相关后的维纳滤波
subplot(224),imshow(fr3,[]);
title('使用自相关函数的维纳滤波后');

fr = deconvreg(g,NOTSEPOWER,RANGE)

%（约束最小二乘方）正则滤波
clc
clear
f=checkerboard(8);%直接生产黑白棋盘，每个小正方形一边的像素点为8个

PSF=fspecial('motion','circular');%具体motion滤波器是什么还不是很清楚。

noise=imnoise(zeros(size(f)),0.001);%产生高斯噪声，且后面要用到该噪声图像

g=gb+noise;
figure,subplot(131),[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');

fr1=deconvreg(g,4);
subplot(132),[]);
title('噪声功率为4的正则滤波后结果');

fr2=deconvreg(g,0.4,[1e-7,1e7]);
subplot(133),[]);
title('带搜索范围的正则滤波后结果');