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

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

clcclearr=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');%结果显示如下：

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

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

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

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

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

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

clcclearC1=[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');%结果如下所示：

%% 估计噪声参数clcclearf=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，总像素点数为npixfigure,bar(p,1);%绘制出条形图title('交互式选取区域产生的直方图');

axis tight

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

%% 掩膜的使用方法clcclearf=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区域大小，列向量%其运算结果为：

%% 空间噪声滤波clcclearf=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;%在原始图像对应椒点的位置赋值为0figure,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;%在原始图像对应盐点的位置赋值为1figure,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('最小值滤波后的结果');

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

g=imnoise(f,'salt & pepper',0.25);%噪声点有白有黑，因为这是用的imnoise,不是imnoise2和imnoise3figure,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的自适应中值滤波器滤波后的图像');

%% 模糊噪声图像建模clcclearf=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('使用自相关函数的维纳滤波后');%其运行结果如下：

%% 约束最小二乘法（正则）滤波clcclearf=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算法的迭代非线性复原clcclearf=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次后非线性复原的图像');%其运行结果如下：

%% 盲目去卷积clcclearPSF=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次后的结果');%其运行结果如下：

%% vistformfwdclcclearTscale=[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);%建立仿射变换tform1figure,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);%建立仿射变换tform2figure,vistformfwd(tform2,[0 100],[0 100]);%其运行结果如下：

%% 图像空间变换clcclearf=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('图像空间变换指定方向和外部填充');%其运行结果如下：

%% 图像配准clcclearg=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('图像配准');%其运行结果如下：