Matlab DIP(瓦)ch4圖像頻域濾波練習

     這一章也是按照網薩雷斯的matlab書練習的，主要講的是圖像的頻域濾波方面的只是，這一次的代碼相對於第3章稍加了些注釋。有幾個函數的功能還不是特別明確。練習代碼如下：

%% fftshift 對數變換,所應用的圖片本身很簡單，就只有黑白2種顏色
clc
clear
f = imread('.\images\dipum_images_ch04\Fig0403(a)(image).tif');
imshow(f)
title('原始圖像')

imfinfo('.\images\dipum_images_ch04\Fig0403(a)(image).tif');%此處如果用Imfinfo(f)就會報錯fft

%沒有居中的傅里葉頻譜
F=fft2(f);%進行二維快速傅里葉變換，其結果和DFT的一樣，只是計算機的計算速度變快了而已，因而叫fft
S=abs(F);%求傅里葉變換后的幅值
figure,subplot(121),imshow(S,[])
title('傅里葉頻譜圖像1');%title函數一定要放在坐標顯示的下一句才有效。

subplot(122),imshow(S),title('傅里葉頻譜圖像2')

%居中的傅里葉頻譜
Fc=fftshift(F);%將頻譜圖像原點移至圖像矩形中間
S1=abs(Fc);
figure,subplot(121),imshow(S1,[]);%加了第二個參數后顯示的圖像正常
subplot(122),imshow(S1);%當沒有第二個參數時，顯示的圖像為豎線加一些孤立的黑點，why？

%使用對數后視覺增強后的傅里葉頻譜
S2=log(1+S1);
imshow(S2,[]);

%%用fftshift對數變換顯示稍復雜的圖片
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0409(a)(bld).tif');
imshow(f);

f=double(f);%其實這句不要試過對后面的變換結果也沒有影響
F=fft2(f);
Fc=fftshift(F);
S=abs(Fc);
S2=log(1+S);%如果沒有這句的話，那么根本看不到細節的圖，所以一定要用對數壓縮增加對比度
figure,imshow(S2,[])

%%試一下ifft功能
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0409(a)(bld).tif');
imshow(f);

f=double(f);
f(1:8,1:8)
F=fft2(f);
f=ifft2(F);%取傅里葉變換后的反傅里葉變換,直接是整數， 並不需要像某些代碼一樣取real部分
f(1:8,1:8)

%%先0擴充再濾波
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0405(a)(square_original).tif');
imshow(f);

[M N]=size(f);
F=fft2(f);%沒有經過0擴充直接計算fft
sig=10;%高斯濾波參數
H=lpfilter('gaussian',M,N,sig);
G=H.*F;%加了.號的乘法表示對應每個元素都相乘，沒加.號的乘法說明是真正的矩陣乘法
g=real(ifft2(G));
figure,imshow(g,[]);

PQ=paddedsize(size(f));%PQ為0擴展的面積，這里設置為與圖像的大小相同？這里size(f)的大小為256*256
Fp=fft2(f,PQ(1),PQ(2));%如果fft2有3個參數的話，則是進行了0擴展了
Hp=lpfilter('gaussian',PQ(1),PQ(2),2*sig);%構造高斯濾波器
Gp=Hp.*Fp;%對0擴展后的圖像進行高斯濾波
gp=real(ifft2(Gp));
figure,imshow(gp,[])'

gpc=gp(1:size(f,1),1:size(f,2));%取gp圖像的1到f的第一維的行，1到f第二維的列部分，即取與圖像大小相同的部分
figure,imshow(gpc,[]);%此時顯示的結果應該與g圖像一樣。

%直接使用空間域濾波
h=fspecial('gaussian',15,7);
gs=imfilter(f,h);
figure,imshow(gs,[])

%%
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0405(a)(square_original).tif');
f=double(f);
imshow(f);

PQ=paddedsize(size(f));%size(f)的大小為256*256
sig=10;
H = lpfilter('gaussian',PQ(1),PQ(2),2*sig); % PQ=[512 512]
figure,mesh(abs(H(1:10:512,1:10:512)));%畫出濾波器的三維空間圖
g=dftfilt(f,H);
figure,imshow(g,[])

%%空間濾波和頻域濾波的比較
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0409(a)(bld).tif');
imshow(f);

F=fft2(f);
S=fftshift(log(1+abs(F)));
S=gscale(S);%用gscale來進行將數據縮放到默認值0~255
figure,imshow(S);

h=fspecial('sobel');%構造sobel空間濾波器
figure,freqz2(h);%查看濾波器的圖形

PQ=paddedsize(size(f));
H=freqz2(h,PQ(1),PQ(2));%將空間濾波器h轉換成頻域濾波器H，但是濾波器的原點在矩陣的中心處
H1=ifftshift(H);%原點遷移到左上角

figure,mesh(abs(H1(1:20:1200,1:20:1200)));
figure,subplot(121),imshow(abs(H),[]);
subplot(122),imshow(abs(H1),[]);

gs=imfilter(double(f),h);%空間濾波，其實就是邊緣檢測
gf=dftfilt(f,H);%頻域濾波
subplot(121),imshow(gs,[]),subplot(122),imshow(gf,[]);%將2種濾波結果顯示出來

figure,imshow(abs(gs)>0.2*abs(max(gs(:))));%只顯示最大值的20%的邊緣
figure,imshow(abs(gf)>0.2*abs(max(gf(:))));

d=abs(gs-gf);
a=max(d(:))
b=min(d(:))

%%構造低通濾波器
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0413(a)(original_test_pattern).tif');
imshow(f)

F1=fft2(f);
imshow(log(1+abs(fftshift(F1))),[]);%直接顯示取對數后的傅里葉變換

PQ=paddedsize(size(f));
[U V]=dftuv(PQ(1),PQ(2));%這個函數還真不明白什么意思，help中解釋的是計算網格頻率矩陣U和V，這2個矩陣是用來頻率濾波的
D0=0.05*PQ(2);

F=fft2(f,PQ(1),PQ(2));%用0擴充大小為PQ(1)*PQ(2)，然后fft變換
figure,imshow(log(1+abs(fftshift(F))),[]);%顯示0擴充后的fft變換圖

H=exp(-(U.^2+V.^2)/(2*(D0^2)));%構造頻域濾波算子
figure,imshow(fftshift(H),[]);%顯示平移后濾波器算子
g=dftfilt(f,H)
figure,imshow(g,[]);

%%繪制一個低通濾波器的mesh圖
clc
clear
H1=lpfilter('gaussian',500,500,50)%構造一個中心點在左上角的高斯低通濾波器,size為500*500
H2=fftshift(H1);%將中心點平移至矩陣中心
mesh(H1(1:10:500,1:10:500)),figure,mesh(H2(1:10:500,1:10:500));%繪出2者的mesh圖
axis([0 50 0 50 0 1]);%xy軸范圍0~50，z軸范圍0~1
figure,subplot(121),imshow(H1,[]),subplot(122),imshow(H2,[]);%繪出2者的灰度圖

%%繪制一個低通濾波器的mesh圖
clc
clear
H=fftshift(hpfilter('ideal',500,500,100));%構造理想的高通濾波器
mesh(H(1:10:500,1:10:500));
axis([1 50 1 50 0 1]);
colormap([0 0 0]);%mesh網格全部用黑色畫
axis off%關掉坐標軸顯示
grid off%關掉網格顯示
figure,imshow(H,[])

%%使用高通濾波器
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0413(a)(original_test_pattern).tif');
imshow(f)
title('原始圖像')

PQ=paddedsize(size(f));
D0=0.05*PQ(1);
H=hpfilter('gaussian',PQ(1),PQ(2),D0);
g=dftfilt(f,H);
figure,imshow(g,[])

%%將高通濾波和直方圖均衡結合起來
clc
clear
f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_original).tif');
imshow(f)
title('原始圖像')

PQ = paddedsize(size(f));
D0 = 0.05*PQ(1);
HBW=hpfilter('btw',PQ(1),PQ(2),D0,2);%構造高通Butterworth濾波器，截斷頻率為D0
H=0.5+2*HBW;
gdw=dftfilt(f,HBW);
gbw=gscale(gdw);%將灰度值自動縮放到0~255之間
figure,subplot(121),imshow(gdw,[]);
title('高通濾波后圖像');

gbe=histeq(gbw,256);%直方圖均衡化
subplot(122),imshow(gbe,[]);
title('高通濾波且直方圖均衡化后圖像');

ghf=dftfilt(f,H);%H是HBW的線性變換
ghf=gscale(ghf);
figure,subplot(121),imshow(ghf,[]);
title('高頻強調濾波后圖像');

ghe=histeq(ghf,256);
subplot(122),imshow(ghe,[]);
title('高頻強調且均衡化后圖像');

%%fftshift和ifftshift的加深理解
clc
clear
A=[2 0 0 1
    0 0 0 0
    0 0 0 0
    3 0 0 4]
B=fftshift(A)%中心點平移
C=fftshift(fftshift(A))%還原成A
D=ifftshift(fftshift(A))%也同樣還原成A
E=ifftshift(A)%平移結果和B一樣

注：

     freqz2 生成的濾波器原點在正中央；lpfilter（低通）生成的濾波器原點在左上角；hpfilter（高通）生成的濾波器原點在左上角

     對0擴充圖像的理解還不是很透，也就是函數paddedsize不是很清楚，返回值的意義也不是很了解；對gscale函數的功能不是很了解，好像是直接將圖像值壓縮到指定的范圍，相信以后會慢慢了解的!

    估計我的注釋過程有的可能有理解錯誤，希望一起交流！