Matlab DIP(瓦)ch5图像复原练习

      这一章内容比较多点，主要将的是图像复原部分，包过线性复原和非线性复原，最好还有一些图像变换方面的练习。这次练习相对上一章把一些比较重要的图片结果贴出来了，希望与大家一起交流！

%% 模拟产生各种噪声
clc
clear
%注意此处函数imnoise2与imnoise不同，imnoise是对一幅图像加噪声
r=imnoise2('gaussian',100000,1,0,1);%imnoise2是产生噪声矩阵，这里是产生高斯噪声，矩阵大小为10000*1
bins=100;                           %均值为0，方差为1
hist(r,bins);%将r矩阵中的数值直方图表示出来，共100个bin
title('guassian');
%结果显示如下：

clc
clear
r=imnoise2('uniform',100000,1,0,1);%同上，此处产生的是0~1的均匀分布噪声，矩阵大小为100000*1
bins=100;
hist(r,bins);
title('uniform');
%结果显示如下：

clc
clear
r=imnoise2('salt & pepper',1000,1,0.1,0.27);%注意这里的salt & pepper中间记得留空格。
bins=100;                                   %为0的概率0.1，为1的概率0.27。为什么加起来不等于1呢？
hist(r,bins);
title('salt & pepper');
%结果显示如下：

clc
clear
r=imnoise2('lognormal',100000,1);%产生对数正态噪声，大小100000*1，偏移值默认为1，形状参数默认0.25
bins=100;
hist(r,bins);
title('lognormal');
%结果显示如下：

clc
clear
r=imnoise2('rayleigh',100000,1,0,1);%rayleigh噪声，该噪声好像是对整体矢量合成的一种分布
bins=100;
hist(r,bins);
title('rayleigh');
%结果显示如下：

clc
clear
r=imnoise2('exponential',100000,1);%指数分布的参数默认为1
bins=100;
hist(r,bins);
title('exponential');
%结果显示如下：

clc
clear
r=imnoise2('erlang',100000,1);%erlang分布，默认值为A=2，B=5。此分布与指数分布和gamma分布有一定的联系
bins=100;
hist(r,bins);
title('erlang');
%结果显示如下：

%%imnoise3练习1
clc
clear
C=[0 64;0 128;32 32;64 0;128 0;-32,32];
[r,R,S]=imnoise3(512,521,C);%这句老是出现错误！！
imshow(S,[]);title('6个指定冲击的正弦周期频谱');
figure,imshow(r,[]);title('6个相应的正弦周期模式');

%%imnoise3练习2
clc
clear
C1=[6 32];
[r,R,S]=imnoise3(512,512,C1);%这个函数功能没怎么弄清楚还！
imshow(S,[]);title('1个指定冲击的正弦噪声周期频谱1');
figure,imshow(r,[]);title('1个相应的正弦噪声周期模式1');
%结果如下所示：

clc
clear
C1=[-2 -2];
[r,R,S]=imnoise3(512,512,C1);%这个函数功能没怎么弄清楚还！
imshow(S,[]);title('1个指定冲击的正弦噪声周期频谱2');
figure,imshow(r,[]);title('1个相应的正弦噪声周期模式2');
%结果如下所示：

clc
clear
C1=[6 32;-2 2];
A=[1 5];
[r,R,S]=imnoise3(512,512,C1,A);%这个函数功能没怎么弄清楚还！
imshow(1-S,[]);title('使用非默认的不同振幅指定冲击的正弦噪声周期频谱1');
figure,imshow(r,[]);title('使用非默认的不同振幅]相应的正弦噪声周期模式2');
%结果如下所示：

%% 估计噪声参数
clc
clear
f=imread('.\images\dipum_images_ch05\Fig0504(a)(noisy_image).tif');
imshow(f),title('原始含噪图像');

[B,c,r]=roipoly(f);%函数roipoly允许用户在图像f中用鼠标标出多边形来，
figure,imshow(B);  %其对应的顶点坐标列与行分别存入向量c和r中，B为其二值图像，里面为1，外面为0

[p,npix]=histroi(f,c,r);%此函数将图像多边形内部分转换成直方图p，总像素点数为npix
figure,bar(p,1);%绘制出条形图
title('交互式选取区域产生的直方图');

axis tight

[v,unv]=statmoments(p,2);%对感兴趣图像区域求均值v和方差unv
X=imnoise2('gaussian',npix,1,unv(1),sqrt(unv(2)));%产生高斯噪声，大小均值方差与感兴趣区域图像的一样
figure,hist(X,150);
axis tight

%% 掩膜的使用方法
clc
clear
f=imread('.\images\dipum_images_ch05\Fig0504(a)(noisy_image).tif');
imshow(f);

[B,c,r]=roipoly(f);%函数roipoly允许用户在图像f中用鼠标标出多边形来
roi=f(B);%将图像f与图像B进行掩膜操作，得到的不规则图像为roi

size_f=size(f)%为原始图像的大小
class_f=class(f)%为原始图像的类型
size_B=size(B)%B图像大小,其实与原始图像大小一样
class_B=class(B)%类似
size_roi=size(roi)%roi区域大小，列向量
%其运算结果为：

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

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

R=imnoise2('salt & pepper',M,N,0,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,3,1.5);%反调和滤波，尺寸大小为3*3，Q默认为1.5,正的Q过滤胡椒噪声
figure,imshow(fp);
title('正Q反调和滤波器滤波的结果')

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

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

fsmin=spfilt(gs,'min',3,3);
figure,imshow(fsmin);
title('最小值滤波后的结果');

%% 自适应中值滤波器
clc
clear
f=imread('.\images\dipum_images_ch05\Fig0507(a)(checkeboard8_pixeldup_8).tif');
imshow(f),title('原始图像');

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

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

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

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

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

noise=imnoise(zeros(size(f)),'gaussian',0,0.001);%产生高斯噪声，且后面要用到该噪声图像
figure,imshow(noise,[]);
title('高斯纯噪声')
%其运行结果如下：

g=gb+noise;
figure,imshow(g,[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');
%其运行结果如下：

%%使用deconvwnr函数复原模糊噪声图像
fr1=deconvwnr(g,PSF);%维纳滤波，去模糊化
figure,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,PSF,R);%信噪比为常数的参数维纳滤波器
figure,imshow(fr2,[]);
title('常数比率信噪比维纳滤波器');
%其运行结果如下：

NCORR=fftshift(real(ifft(Sn)));%噪声的自相关函数，why？
ICORR=fftshift(real(ifft(Sf)));%信号的自相关函数，why？
fr3=deconvwnr(g,PSF,NCORR,ICORR);%自相关后的维纳滤波
figure,imshow(fr3,[]);
title('使用自相关函数的维纳滤波后');
%其运行结果如下：

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

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

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

g=gb+noise;
imshow(g,[]);%模糊加噪声后的图像
title('模糊加噪声后的图像');
%其运行结果如下：

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

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

%% 使用L-R算法的迭代非线性复原
clc
clear
f=checkerboard(8);
imshow(pixeldup(f,8));%pixeldup函数是将图像扩大m*n倍，通过复制每个像素点m*n次。
title('原始图像');
%其运行结果如下：

PSF=fspecial('gaussian',7,10);%为什么这个时候的PSFsize是1*3呢，按理说应该是7*7的
SD=0.01;
g=imnoise(imfilter(f,PSF),'gaussian',0,SD^2);%加了2次高斯模糊
figure,imshow(pixeldup(g,8));
title('模糊加噪的图像');

%其运行结果如下：

DAMPAR=10*SD;

LIM=ceil(size(PSF,1)/2);%ceil函数是求最接近的整数
WEIGHT=zeros(size(g));
WEIGHT(LIM+1:end-LIM,LIM+1:end-LIM)=1;

NUMIT=5;
f5=deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);%采用LR算法非线性复原
figure,imshow(pixeldup(f5,8));
title('使用LR算法迭代5次后非线性复原的图像');
%其运行结果如下：

NUMIT=10;
f10=deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);%采用LR算法非线性复原
figure,imshow(pixeldup(f10,8));
title('使用LR算法迭代10次后非线性复原的图像');
%其运行结果如下：

NUMIT=20;
f20=deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);%采用LR算法非线性复原
figure,imshow(pixeldup(f20,8));
title('使用LR算法迭代20次后非线性复原的图像');
%其运行结果如下：

NUMIT=100;
f100=deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);%采用LR算法非线性复原
figure,imshow(pixeldup(f100,8));
title('使用LR算法迭代100次后非线性复原的图像');
%其运行结果如下：

NUMIT=1000;
f1000=deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT);%采用LR算法非线性复原
figure,imshow(pixeldup(f1000,8));
title('使用LR算法迭代1000次后非线性复原的图像');
%其运行结果如下：

%% 盲目去卷积
clc
clear
PSF=fspecial('gaussian',7,10);
imshow(pixeldup(PSF,73),[]);
title('原始图像');
%其运行结果如下：

f=checkerboard(8);
SD=0.01;
g=imnoise(imfilter(f,PSF),'gaussian',0,SD^2);

INITPSF=ones(size(PSF));%PSF的初始估计矩阵
DAMPAR=10*SD;
LIM=ceil(size(PSF,1)/2);
WEIGHT=zeros(size(g));
WEIGHT(LIM+1:end-LIM,LIM+1:end-LIM)=1;

NUMIT=5;
[fr,PSFe]=deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);%使用盲去卷积迭代复原，返回的PSFe是点扩散函数
figure,imshow(pixeldup(PSFe,73),[]);
title('使用盲去卷积估计PSF迭代5次后的结果');

%其运行结果如下：

NUMIT=10;
[fr,PSFe]=deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);%使用盲去卷积迭代复原，返回的PSFe是点扩散函数
figure,imshow(pixeldup(PSFe,73),[]);
title('使用盲去卷积估计PSF迭代10次后的结果');
%其运行结果如下：

NUMIT=20;
[fr,PSFe]=deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);%使用盲去卷积迭代复原，返回的PSFe是点扩散函数
figure,imshow(pixeldup(PSFe,73),[]);
title('使用盲去卷积估计PSF迭代20次后的结果');
%其运行结果如下：

NUMIT=50;
[fr,PSFe]=deconvblind(g,INITPSF,NUMIT,DAMPAR,WEIGHT);%使用盲去卷积迭代复原，返回的PSFe是点扩散函数
figure,imshow(pixeldup(PSFe,73),[]);
title('使用盲去卷积估计PSF迭代50次后的结果');
%其运行结果如下：

%% vistformfwd
clc
clear
Tscale=[1.5 0 0;0 2 0;0 0 1];%仿射变换矩阵尺度部分
Trotation=[cos(pi/4) sin(pi/4) 0%仿射变换矩阵旋转部分
            -sin(pi/4) cos(pi/4) 0
            0 0 1];
T1=Tscale*Trotation;%仿射变换矩阵
tform1=maketform('affine',T1);%建立仿射变换tform1
figure,vistformfwd(tform1,[0 100],[0 100]);
%其运行结果如下：

Tshear=[1 0 0;0.2 1 0;0 0 1];%仿射矩阵的水平剪枝部分
T2=Tscale*Trotation*Tshear;
tform2=maketform('affine',T2);%建立仿射变换tform2
figure,vistformfwd(tform2,[0 100],[0 100]);
%其运行结果如下：

%% 图像空间变换
clc
clear
f=checkerboard(50);%8*8的checkboard，每个小正方形有50个像素
imshow(f);
title('图像空间变换原始图');
%其运行结果如下：

s=0.8;
theta=pi/6;
T=[s*cos(theta) s*sin(theta) 0 
    -s*sin(theta) s*cos(theta) 0 
    0 0 1];
tform=maketform('affine',T);
g=imtransform(f,tform);
figure,imshow(g,[]);
title('图像空间变换1');
%其运行结果如下：

g2=imtransform(f,tform,'nearest');
figure,imshow(g2,[]);
title('图像空间变换最近邻插值');
%其运行结果如下：

g3=imtransform(f,tform,'FillValue',0.5);
figure,imshow(g3,[]);
title('图像空间变换外部0.5填充');
%其运行结果如下：

T2 = [1 0 0; 0 1 0; 50 50 1];
tform2 = maketform('affine',T2);
g4 = imtransform(f,tform2);
figure,imshow(g4,[]);
title('图像空间变换平移');
%其运行结果如下：

g5 = imtransform(f,tform2,'XData',[1 500],'YData',[1 500],...
                 'FillValue',0.5);
figure,imshow(g5,[]);
title('图像空间变换指定方向和外部填充');
%其运行结果如下：

%% 图像配准
clc
clear
g=imread('.\images\dipum_images_ch05\Fig0515(a)(base-with-control-points).tif');
imshow(g,[]);
title('原始图像');
%其运行结果如下：

basepoints=[83 81;450 56;43 293;249 392;436 442];
inputpoints=[68 66;375 47;42 286;275 434;523 532];
tform=cp2tform(inputpoints,basepoints,'projective');
gp=imtransform(g,tform,'XData',[1 502],'YData',[1 502]);
figure,imshow(gp,[]);
title('图像配准');
%其运行结果如下：

 

      从上面的练习过程可以看出，在进行迭代复原的过程中，并不是迭代次数越多就越明显。另外如果我们已经对噪声和图像谱的知识有足够的了解的前提下。Wiener滤波结果要好得多。如果没有这些信息，则用“约束的最小二乘法（正则）”滤波器和Wiener滤波基本差不多。

      欢迎交流！