這一章內容比較多點，主要將的是圖像復原部分，包過線性復原和非線性復原，最好還有一些圖像變換方面的練習。這次練習相對上一章把一些比較重要的圖片結果貼出來了，希望與大家一起交流！
%% 模擬產生各種噪聲
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濾波基本差不多。
歡迎交流！