圖像恢複及濾波處理
1.基本概念：
a.圖像恢複是通過計算機處理，對品質下降的圖像加以重建或恢複的處理過程。因錄影機與物體相對運動、系統誤差、畸變、噪聲等因素的影響，使圖像往往不是真實景物的完善映像。在圖像恢複中，需建立造成圖像品質下降的退化模型，然後運用相反過程來恢複原來圖像，并運用一定準則來判定是否得到圖像的最佳恢複。在遙感圖像進行中，為消除遙感圖像的失真、畸變，恢複目标的反射波譜特性和正确的幾何位置，通常需要對圖像進行恢複處理，包括輻射校正、大氣校正、條帶噪聲消除、幾何校正等内容。
b.圖像恢複采用的方法：
1）空間域濾波：均值濾波，順序統計濾波； 2）頻率域濾波：帶通、帶阻濾波，陷波濾波。
濾波器被廣泛地用于圖象的預處理，抑制圖象噪聲,增強對比度，以及強化圖象的邊沿特征。運用較為廣泛的線性濾波器如平均值濾波器，能較好地抑制圖象中的加性噪聲。 但是，線性濾波器會引起圖象的鈍化或模糊，使得圖象中物體邊界産生位移。特别是，在圖象受到乘性噪聲或脈沖噪聲的幹擾，如超音波及雷達成像中普遍存在的斑點噪聲，線性濾波器就不能取得預期的效果。中值濾波器，就像其名字一樣，是用該像素的相鄰像素的灰階中值來代替該像素的值，是一種非線性濾波器。
c.椒鹽噪聲：噪聲脈沖可正可負，正脈沖是以白點（鹽點）出現，負脈沖是以黑點（胡椒點）出現。
d.算術均值濾波器和幾何均值濾波器（尤其是後者）更适合處理高斯或均勻等随機噪聲。
諧波均值濾波器：對于“鹽”噪聲效果好，不适用“胡椒”噪聲，善于處理高斯噪聲那樣的脈沖噪聲。
缺點是：必須知道噪聲是暗噪聲還是亮噪聲，以便選擇合适的Q号。
逆諧波均值濾波器：Q為其階數，當Q為正數，用于消除“胡椒”噪聲；當Q為負數，用于消除“鹽”噪聲；
當Q=0退變為算術均值濾波器，當Q=-1時，退變為諧波均值濾波器。
e.順序統計濾波器：
中值濾波器:處理脈沖型加性噪聲，重複使用中值濾波器處理圖像會使圖像模糊，盡可能保持希望的處理次數。
最大值：除去圖像中的“胡椒”噪聲，但是會從黑色物體的邊緣移走一些黑色像素。
最大值：除去圖像中的“鹽”噪聲，但是會從白色物體的邊緣移走一些白色像素。
修正的阿爾法均值濾波器：d值可以取0到mn-1之間的任意數，當d=0時，退變為算術均值濾波器,
當d=（mn-1）/2退變為中值濾波器，d取其他值時，在包含多種噪聲的情況下非常适用。
f.巴特沃斯IIR濾波器的設計:
在MATLAB下，設計巴特沃斯IIR濾波器可使用butter函數。
Butter函數可設計低通、高通、帶通和帶阻的數字和模拟IIR濾波器，其特性為使通帶内的幅度響應最大限度地平坦，但同時損失截止頻率處的下降斜度。在期望通帶平滑的情況下，可使用butter函數。
butter函數的用法為：[b,a]=butter(n,Wn,/ftype/)
其中n代表濾波器階數，Wn代表濾波器的截止頻率，這兩個參數可使用buttord函數來确定。buttord函數可在給定濾波器性能的情況下，求出巴特沃斯濾波器的最小階數n，同時給出對應的截止頻率Wn。
buttord函數的用法為：[n,Wn]= buttord(Wp,Ws,Rp,Rs)
其中Wp和Ws分别是通帶和阻帶的拐角頻率（截止頻率），其取值範圍為0至1之間。當其值為1時代表采樣頻率的一半。Rp和Rs分别是通帶和阻帶區的波紋系數。
g.不同類型（高通、低通、帶通和帶阻）濾波器對應的Wp和Ws值遵循以下規則：
(1) 高通濾波器：Wp和Ws為一進制矢量且Wp>Ws；
(2) 低通濾波器：Wp和Ws為一進制矢量且Wp<Ws；
(3) 帶通濾波器：Wp和Ws為二進制矢量且Wp<Ws，如Wp=[0.2,0.7],Ws=[0.1,0.8];
(4) 帶阻濾波器：Wp和Ws為二進制矢量且Wp>Ws，如Wp=[0.1,0.8],Ws=[0.2,0.7]。
h.契比雪夫I型IIR濾波器的設計
在期望通帶下降斜率大的場合，應使用橢圓濾波器或契比雪夫濾波器。在MATLAB下可使用cheby1函數設計出契比雪夫I型IIR濾波器。
cheby1函數可設計低通、高通、帶通和帶阻契比雪夫I型濾IIR波器，其通帶内為等波紋，阻帶内為單調。契比雪夫I型的下降斜度比II型大，但其代價是通帶内波紋較大。
cheby1函數的用法為：[b,a]=cheby1(n,Rp,Wn,/ftype/)
在使用cheby1函數設計IIR濾波器之前，可使用cheblord函數求出濾波器階數n和截止頻率Wn。cheblord函數可在給定濾波器性能的情況下，選擇契比雪夫I型濾波器的最小階和截止頻率Wn。
cheblord函數的用法為：[n,Wn]=cheblord(Wp,Ws,Rp,Rs)
其中Wp和Ws分别是通帶和阻帶的拐角頻率（截止頻率），其取值範圍為0至1之間。當其值為1時代表采樣頻率的一半。Rp和Rs分别是通帶和阻帶區的波紋系數。
2.常用函數：
Spfilt：空間域濾波；
Lpfilter：低通濾波；
Imnoise2：生成指定分布模型的噪聲；
Imnoise3：生成周期噪聲 。
3.應用：
（1）處理過程：
• 無padding過程
直接在頻率域産生于圖像大小一緻的濾波器。
• 有padding過程
頻域濾波函數乘以
的反變換，用0延拓，然後再進行傅裡葉正變換。
a.無padding過程:
顯示原圖像：
f=imread(windmill_noise.png');
figure,subplot(2,2,1),imshow(f);
title(‘原圖像');
F=fft2(f);
[M N]=size(f);
f=double(f);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,2),imshow(s,[]);
title(‘原圖像頻譜');<span style="color:#cc66cc;">
理想帶阻濾波器濾波效果：
[u1,v1]=dftuv(M,N);
d1=sqrt(u1.^2+v1.^2);
w=5;
D0=15;
ba= d1<(D0-w/2) | d1>(D0+w/2);
H=double(ba);
g=real(ifft2(H.*fft2(f)));
subplot(2,2,3),imshow(g,[]);
title(‘理想帶阻');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]);
title(理想帶阻頻譜');
巴特沃斯帶阻濾波：
n=3；
H = 1./(1 + (w*d1./(d1.^2-D0^2)).^(2*n));
g=real(ifft2(H.*fft2(f)));
figure,subplot(2,2,1),imshow(g,[]);
title(‘btr');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,2),imshow(s,[]);
title(‘btr頻譜圖');
高斯帶阻濾波：
H = 1-exp(-1/2*(((d1.^2)-D0^2)./(d1*w)).^2);
g=real(ifft2(H.*fft2(f)));
subplot(2,2,3),imshow(g,[]);
title(‘高斯帶阻');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]);
title(‘高斯帶阻頻譜');
（2）主要MATLAB源碼：
① dftfilt.m
function g=dftfilt(f,H)
F=fft2(f,size(H,1),size(H,2));
g=real(ifft2(H.*F));
g=g(1:size(f,1),1:size(f,2));
② dftuv.m
function [u,v]=dftuv(m,n)
u=0:(m-1);
v=0:(n-1);
idx=find(u>m/2);
u(idx)=u(idx)-m;
idy=find(v>n/2);
v(idy)=v(idy)-n;
[v,u]=meshgrid(v,u);
③ pad1.m
%本檔案包含paddedsize.m和dftuv.m及dftfilt.m三個檔案
f=imread('D:\test\windmill_noise.png');
figure,subplot(2,2,1),imshow(f); title('原圖像');
F=fft2(f);
[M N]=size(f);
f=double(f);
%fc=ifft2(fft2(f));
fc=fftshift(F);
s=log(1+abs(fc));
% subplot(2,2,2),imshow(s,[]); title('原圖像頻譜'); %原圖像頻譜
%##########################
PQ=paddedsize(size(f));
[u,v]=dftuv(PQ(1),PQ(2));
[u1,v1]=dftuv(M,N); %沒有padding
d1=sqrt(u1.^2+v1.^2);
w=10;%帶寬 理想帶阻濾波
D0=15;
ba= d1<(D0-w/2) | d1>(D0+w/2);
H=double(ba);
subplot(2,2,2),imshow(fftshift(H),[]); title('H');
%H = 1-exp(-1/2*(((d1.^2)-D0^2)./(d1*w)).^2);%高斯帶阻濾波器
% H=fftshift(H);
% H=ifftshift(H);
% h=(-1).^(u1+v1);
% hh=ifft2(h);
H=real(ifft2(H));
% subplot(2,2,2),imshow(H,[]); title('H');
H1=imag(ifft2(H)); H=H+1i*H1;
H=ifftshift(H);
% H=H.*hh;
H=fft2(H,900,1200);
% subplot(2,2,2),imshow(ifft(real(H)),[]); title('H');
F=fft2(f,900,1200);
% F=fftshift(F);
% H=H(1:size(f,1),1:size(f,2));
% F=F(1:size(f,1),1:size(f,2));
tmp=H.*F;
g=real(ifft2(H.*F));
%subplot(2,2,3),imshow(ifftshift(log(1+abs(F.*H))),[]); title('無padding高斯帶阻圖像頻譜');
%g=g(226:675,301:900);
% g=g(1:size(f,1),1:size(f,2));
subplot(2,2,3),imshow(g,[]); title('無padding高斯帶阻');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]); title('無padding高斯帶阻圖像頻譜');
④ pad2.m
%本檔案包含paddedsize.m和dftuv.m及dftfilt.m三個檔案
f=imread('D:\test\windmill_noise.png');
figure,subplot(2,2,1),imshow(f); title('原圖像');
F=fft2(f);
[M N]=size(f);
f=double(f);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,2),imshow(s,[]); title('原圖像頻譜'); %原圖像頻譜
%##########################
% PQ=paddedsize(size(f));
% [u,v]=dftuv(PQ(1),PQ(2));
[u1,v1]=dftuv(M,N); %沒有padding
d1=sqrt(u1.^2+v1.^2);
%d0=0.05*PQ(2);
%d0=0.1*N;
d0=15;
% F=fft2(f,PQ(1),PQ(2));
% H=1-(1-exp(-(u1.^2+v1.^2)/(2*(d0^2)))); %高斯高通濾波
% %g=dftfilt(f,H);
% g=real(ifft2(H.*fft2(f)));
% subplot(2,2,3),imshow(g,[]); title('高斯high通濾波');
% F=fft2(g);
% fc=fftshift(F);
% s=log(1+abs(fc));
% subplot(2,2,4),imshow(s,[]); title('高斯濾波後的圖像頻譜'); %高斯濾波後的圖像頻譜
% %###########################
% d=sqrt(u.^2+v.^2);
% H=1-1./(1+(d1/d0).^2); %巴特沃茲高通濾波
% %g=dftfilt(f,H);
% g=real(ifft2(H.*fft2(f)));
% figure,subplot(2,2,1),imshow(g,[]); title('巴特沃茲high通濾波');
% F=fft2(g);
% fc=fftshift(F);
% s=log(1+abs(fc));
% subplot(2,2,2),imshow(s,[]); title('巴特沃茲high通濾波後的圖像頻譜'); %巴特沃茲低通濾波後的圖像頻譜
% %##########################
% h=double(d1<=d0); %理想高通濾波
% %g=dftfilt(f,H);
% g=real(ifft2(H.*fft2(f)));
% H=1-double(d1<=d0);
% subplot(2,2,3),imshow(g,[]); title('理想high通濾波');
% F=fft2(g);
% fc=fftshift(F);
% s=log(1+abs(fc));
% subplot(2,2,4),imshow(s,[]); title('理想high通濾波後的圖像頻譜'); %理想低通濾波後的圖像頻譜
%#####################
w=5;%帶寬 理想帶阻濾波
% % D0=sqrt(209*209+298*298);
D0=15;
ba= d1<(D0-w/2) | d1>(D0+w/2);
H=double(ba);
%g=dftfilt(f,H);
g=real(ifft2(H.*fft2(f)));
subplot(2,2,3),imshow(g,[]); title('理想帶阻濾波');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]); title('理想帶阻濾波後的圖像頻譜'); %理想帶阻濾波後的圖像頻譜
%###########################
n=3; %3階巴特沃斯濾波器
H = 1./(1 + (w*d1./(d1.^2-D0^2)).^(2*n));
%g=dftfilt(f,H);
g=real(ifft2(H.*fft2(f)));
figure,subplot(2,2,1),imshow(g,[]); title('巴特沃斯帶阻濾波');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,2),imshow(s,[]); title('巴特沃斯帶阻濾波圖像頻譜'); %巴特沃斯帶阻濾波後的圖像頻譜
%#########################
H = 1-exp(-1/2*(((d1.^2)-D0^2)./(d1*w)).^2);%高斯帶阻濾波器
%g=dftfilt(f,H);
g=real(ifft2(H.*fft2(f)));
subplot(2,2,3),imshow(g,[]); title('高斯帶阻濾波器');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]); title('高斯帶阻濾波後的圖像頻譜'); %高斯帶阻濾波後的圖像頻譜
% %###########################
% v0=15;
% u0=0;
% d1=((u-M/2-u0).^2 + (v-N/2-v0).^2).^(1/2);
% d2=((u-M/2-u0).^2 + (v-N/2+v0).^2).^(1/2);
% H=double(d0<d1 & d0<d2);%理想陷波帶阻濾波器
% g=dftfilt(f,H);
% subplot(2,2,3),imshow(g,[]);
% F=fft2(g);
% fc=fftshift(F);
% s=log(1+abs(fc));
% subplot(2,2,4),imshow(s,[]); %理想陷波帶阻濾波後的圖像頻譜
% %H=1/(1+(D0^2/(d1.*d2)).^2);
⑤ pad3.m
%本檔案包含paddedsize.m和dftuv.m及dftfilt.m三個檔案
f=imread('D:\test\windmill_noise.png');
figure,subplot(2,2,1),imshow(f); title('原圖像');
F=fft2(f);
[M N]=size(f);
f=double(f);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,2),imshow(s,[]); title('原圖像頻譜'); %原圖像頻譜
%##########################
[u1,v1]=dftuv(M,N); %沒有padding
d1=sqrt(u1.^2+v1.^2);
d0=15;
%#####################
w=5;%帶寬 理想帶阻濾波
D0=15;
ba= d1<(D0-w/2) | d1>(D0+w/2);
H=double(ba);
g=real(ifft2(H.*fft2(f)));
subplot(2,2,3),imshow(g,[]); title('理想帶阻濾波');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]); title('理想帶阻濾波後的圖像頻譜'); %理想帶阻濾波後的圖像頻譜
%###########################
n=3; %3階巴特沃斯濾波器
H = 1./(1 + (w*d1./(d1.^2-D0^2)).^(2*n));
g=real(ifft2(H.*fft2(f)));
figure,subplot(2,2,1),imshow(g,[]); title('巴特沃斯帶阻濾波');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,2),imshow(s,[]); title('巴特沃斯帶阻濾波圖像頻譜'); %巴特沃斯帶阻濾波後的圖像頻譜
%#########################
H = 1-exp(-1/2*(((d1.^2)-D0^2)./(d1*w)).^2);%高斯帶阻濾波器
g=real(ifft2(H.*fft2(f)));
subplot(2,2,3),imshow(g,[]); title('高斯帶阻濾波器');
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]); title('高斯帶阻濾波後的圖像頻譜'); %高斯帶阻濾波後的圖像頻譜
⑥ paddedsize.m
function PQ = paddedsize(AB, CD, PARAM)
%PADDEDSIZE Computes padded sizes useful for FFT-based filtering.
% PQ = PADDEDSIZE(AB), where AB is a two-element size vector,
% computes the two-element size vector PQ = 2*AB.
%
% PQ = PADDEDSIZE(AB, 'PWR2') computes the vector PQ such that
% PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX(AB).
%
% PQ = PADDEDSIZE(AB, CD), where AB and CD are two-element size
% vectors, computes the two-element size vector PQ. The elements
% of PQ are the smallest even integers greater than or equal to
% AB + CD - 1.
%
% PQ = PADDEDSIZE(AB, CD, 'PWR2') computes the vector PQ such that
% PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX([AB CD]).
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.5 $ $Date: 2003/08/25 14:28:22 $
if nargin == 1
PQ = 2*AB;
elseif nargin == 2 & ~ischar(CD)
PQ = AB + CD - 1;
PQ = 2 * ceil(PQ / 2);
elseif nargin == 2
m = max(AB); % Maximum dimension.
% Find power-of-2 at least twice m.
P = 2^nextpow2(2*m);
PQ = [P, P];
elseif nargin == 3
m = max([AB CD]); % Maximum dimension.
P = 2^nextpow2(2*m);
PQ = [P, P];
else
error('Wrong number of inputs.')
end
以上代碼的處理結果：
4.附錄：
可能用到的其他函數代碼： ① 生成指定分布模型的噪聲: imnoise2.m
function R = imnoise2(type, M, N, a, b)
%IMNOISE2 Generates an array of random numbers with specified PDF.
% R = IMNOISE2(TYPE, M, N, A, B) generates an array, R, of size
% M-by-N, whose elements are random numbers of the specified TYPE
% with parameters A and B. If only TYPE is included in the
% input argument list, a single random number of the specified
% TYPE and default parameters shown below is generated. If only
% TYPE, M, and N are provided, the default parameters shown below
% are used. If M = N = 1, IMNOISE2 generates a single random
% number of the specified TYPE and parameters A and B.
%
% Valid values for TYPE and parameters A and B are:
%
% 'uniform' Uniform random numbers in the interval (A, B).
% The default values are (0, 1).
% 'gaussian' Gaussian random numbers with mean A and standard
% deviation B. The default values are A = 0, B = 1.
% 'salt & pepper' Salt and pepper numbers of amplitude 0 with
% probability Pa = A, and amplitude 1 with
% probability Pb = B. The default values are Pa =
% Pb = A = B = 0.05. Note that the noise has
% values 0 (with probability Pa = A) and 1 (with
% probability Pb = B), so scaling is necessary if
% values other than 0 and 1 are required. The noise
% matrix R is assigned three values. If R(x, y) =
% 0, the noise at (x, y) is pepper (black). If
% R(x, y) = 1, the noise at (x, y) is salt
% (white). If R(x, y) = 0.5, there is no noise
% assigned to coordinates (x, y).
% 'lognormal' Lognormal numbers with offset A and shape
% parameter B. The defaults are A = 1 and B =
% 0.25.
% 'rayleigh' Rayleigh noise with parameters A and B. The
% default values are A = 0 and B = 1.
% 'exponential' Exponential random numbers with parameter A. The
% default is A = 1.
% 'erlang' Erlang (gamma) random numbers with parameters A
% and B. B must be a positive integer. The
% defaults are A = 2 and B = 5. Erlang random
% numbers are approximated as the sum of B
% exponential random numbers.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.5 $ $Date: 2003/10/12 23:37:29 $
% Set default values.
if nargin == 1
a = 0; b = 1;
M = 1; N = 1;
elseif nargin == 3
a = 0; b = 1;
end
% Begin processing. Use lower(type) to protect against input being
% capitalized.
switch lower(type)
case 'uniform'
R = a + (b - a)*rand(M, N);
case 'gaussian'
R = a + b*randn(M, N);
case 'salt & pepper'
if nargin <= 3
a = 0.05; b = 0.05;
end
% Check to make sure that Pa + Pb is not > 1.
if (a + b) > 1
error('The sum Pa + Pb must not exceed 1.')
end
R(1:M, 1:N) = 0.5;
% Generate an M-by-N array of uniformly-distributed random numbers
% in the range (0, 1). Then, Pa*(M*N) of them will have values <=
% a. The coordinates of these points we call 0 (pepper
% noise). Similarly, Pb*(M*N) points will have values in the range
% > a & <= (a + b). These we call 1 (salt noise).
X = rand(M, N);
c = find(X <= a);
R(c) = 0;
u = a + b;
c = find(X > a & X <= u);
R(c) = 1;
case 'lognormal'
if nargin <= 3
a = 1; b = 0.25;
end
R = a*exp(b*randn(M, N));
case 'rayleigh'
R = a + (-b*log(1 - rand(M, N))).^0.5;
case 'exponential'
if nargin <= 3
a = 1;
end
if a <= 0
error('Parameter a must be positive for exponential type.')
end
k = -1/a;
R = k*log(1 - rand(M, N));
case 'erlang'
if nargin <= 3
a = 2; b = 5;
end
if (b ~= round(b) | b <= 0)
error('Param b must be a positive integer for Erlang.')
end
k = -1/a;
R = zeros(M, N);
for j = 1:b
R = R + k*log(1 - rand(M, N));
end
otherwise
error('Unknown distribution type.')
end<span style="color:#3333ff;">
</span>
② 生成周期噪聲: imnoise3.m
function [r, R, S] = imnoise3(M, N, C, A, B)
%IMNOISE3 Generates periodic noise.
% [r, R, S] = IMNOISE3(M, N, C, A, B), generates a spatial
% sinusoidal noise pattern, r, of size M-by-N, its Fourier
% transform, R, and spectrum, S. The remaining parameters are:
%
% C is a K-by-2 matrix with K pairs of freq. domain coordinates (u,
% v) that define the locations of impulses in the freq. domain. The
% locations are with respect to the frequency rectangle center at
% (floor(M/2) + 1, floor(N/2) + 1). The impulse locations are spe-
% cified as increments with respect to the center. For ex, if M =
% N = 512, then the center is at (257, 257). To specify an impulse
% at (280, 300) we specify the pair (23, 43); i.e., 257 + 23 = 280,
% and 257 + 43 = 300. Only one pair of coordinates is required for
% each impulse. The conjugate pairs are generated automatically.
%
% A is a 1-by-K vector that contains the amplitude of each of the
% K impulse pairs. If A is not included in the argument, the
% default used is A = ONES(1, K). B is then automatically set to
% its default values (see next paragraph). The value specified
% for A(j) is associated with the coordinates in C(j, 1:2).
%
% B is a K-by-2 matrix containing the Bx and By phase components
% for each impulse pair. The default values for B are B(1:K, 1:2)
% = 0.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.5 $ $Date: 2004/11/04 22:32:42 $
% Process input parameters.
[K, n] = size(C);
if nargin == 3
A(1:K) = 1.0;
B(1:K, 1:2) = 0;
elseif nargin == 4
B(1:K, 1:2) = 0;
end
% Generate R.
R = zeros(M, N);
for j = 1:K
u1 = M/2 + 1 + C(j, 1); v1 = N/2 + 1 + C(j, 2);
R(u1, v1) = i * (A(j)/2) * exp(i*2*pi*C(j, 1) * B(j, 1)/M);
% Complex conjugate.
u2 = M/2 + 1 - C(j, 1); v2 = N/2 + 1 - C(j, 2);
R(u2, v2) = -i * (A(j)/2) * exp(i*2*pi*C(j, 2) * B(j, 2)/N);
end
% Compute spectrum and spatial sinusoidal pattern.
S = abs(R);
r = real(ifft2(ifftshift(R)));
③ 低通濾波: lpfilter.m
function H = lpfilter(type, M, N, D0, n)
%LPFILTER Computes frequency domain lowpass filters.
% H = LPFILTER(TYPE, M, N, D0, n) creates the transfer function of
% a lowpass filter, H, of the specified TYPE and size (M-by-N). To
% view the filter as an image or mesh plot, it should be centered
% using H = fftshift(H).
%
% Valid values for TYPE, D0, and n are:
%
% 'ideal' Ideal lowpass filter with cutoff frequency D0. n need
% not be supplied. D0 must be positive.
%
% 'btw' Butterworth lowpass filter of order n, and cutoff
% D0. The default value for n is 1.0. D0 must be
% positive.
%
% 'gaussian' Gaussian lowpass filter with cutoff (standard
% deviation) D0. n need not be supplied. D0 must be
% positive.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.8 $ $Date: 2004/11/04 22:33:16 $
% Use function dftuv to set up the meshgrid arrays needed for
% computing the required distances.
[U, V] = dftuv(M, N);
% Compute the distances D(U, V).
D = sqrt(U.^2 + V.^2);
% Begin filter computations.
switch type
case 'ideal'
H = double(D <= D0);
case 'btw'
if nargin == 4
n = 1;
end
H = 1./(1 + (D./D0).^(2*n));
case 'gaussian'
H = exp(-(D.^2)./(2*(D0^2)));
otherwise
error('Unknown filter type.')
end
④ 帶通濾波：bpfilter.m
function H=bpfilter(type, M, N, D0,n,W)
if nargin == 4
n = 1; % Default value of n.
end
% Generate highpass filter.
Hlp = bsfilter(type, M, N, D0, n,W);
H = 1 - Hlp;
⑤ 帶阻濾波：bsfilter.m
function H=bsfilter(type, M, N, D0,n,W)
[U, V] = dftuv(M, N);
% Compute the distances D(U, V).
D = sqrt(U.^2 + V.^2);
% Begin filter computations.
switch type
case 'ideal'
if D < (D0-W/2)|D > (D0+W/2)
H=1;
else
H=0;
end
case 'btw'
if nargin == 4
n = 1;
end
H = 1./(1 + (D./(D.^2-D0^2)).^(2*n));
case 'gaussian'
H = 1-exp-1/2*(((D.^2)-D0^2)./(D*W)).^2;
otherwise
error('Unknown filter type.')
end
⑥ 高通濾波：hpfilter.m
function H = hpfilter(type, M, N, D0, n)
%HPFILTER Computes frequency domain highpass filters.
% H = HPFILTER(TYPE, M, N, D0, n) creates the transfer function of
% a highpass filter, H, of the specified TYPE and size (M-by-N).
% Valid values for TYPE, D0, and n are:
%
% 'ideal' Ideal highpass filter with cutoff frequency D0. n
% need not be supplied. D0 must be positive.
%
% 'btw' Butterworth highpass filter of order n, and cutoff
% D0. The default value for n is 1.0. D0 must be
% positive.
%
% 'gaussian' Gaussian highpass filter with cutoff (standard
% deviation) D0. n need not be supplied. D0 must be
% positive.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.4 $ $Date: 2003/08/25 14:28:22 $
% The transfer function Hhp of a highpass filter is 1 - Hlp,
% where Hlp is the transfer function of the corresponding lowpass
% filter. Thus, we can use function lpfilter to generate highpass
% filters.
if nargin == 4
n = 1; % Default value of n.
end
% Generate highpass filter.
Hlp = lpfilter(type, M, N, D0, n);
H = 1 - Hlp;
⑦ adpmedian.m
function f=adpmedian(g,smax)
if(smax<=1)||(smax/2==round(smax/2))||(smax~=round(smax))
error('smax must be an odd inger>1')
end
[M,N]=size(g);
f=g;
f(:)=0;
alreadyprocessed=false(size(g));
for k=3:2:smax
zmin=ordfilt2(g,1,ones(k,k),'symmetric');
zmax=ordfilt2(g,k*k,ones(k,k),'symmetric');
zmed=medfilt2(g,[k,k],'symmetric');
processusinglevelB=(zmed>zmin)&(zmax>zmed)&~alreadyprocessed;
zB=(g>zmin)&(zmax>g);
outputzxy=processusinglevelB &zB;
outputzmed=processusinglevelB &~zB;
f(outputzxy)=g(outputzxy);
f(outputzmed)=g(outputzmed);
alreadyprocessed=alreadyprocessed|processusinglevelB;
if all( alreadyprocessed(:))
break;
end
end
⑧ 改變圖像的存儲類：changeclass.m
function image = changeclass(class, varargin)
%CHANGECLASS changes the storage class of an image.
% I2 = CHANGECLASS(CLASS, I);
% RGB2 = CHANGECLASS(CLASS, RGB);
% BW2 = CHANGECLASS(CLASS, BW);
% X2 = CHANGECLASS(CLASS, X, 'indexed');
% Copyright 1993-2002 The MathWorks, Inc. Used with permission.
% $Revision: 1.2 $ $Date: 2003/02/19 22:09:58 $
switch class
case 'uint8'
image = im2uint8(varargin{:});
case 'uint16'
image = im2uint16(varargin{:});
case 'double'
image = im2double(varargin{:});
otherwise
error('Unsupported IPT data class.');
end
⑨ 空間域濾波: spfilt.m
function f = spfilt(g, type, m, n, parameter)
%SPFILT Performs linear and nonlinear spatial filtering.
% F = SPFILT(G, TYPE, M, N, PARAMETER) performs spatial filtering
% of image G using a TYPE filter of size M-by-N. Valid calls to
% SPFILT are as follows:
%
% F = SPFILT(G, 'amean', M, N) Arithmetic mean filtering.
% F = SPFILT(G, 'gmean', M, N) Geometric mean filtering.
% F = SPFILT(G, 'hmean', M, N) Harmonic mean filtering.
% F = SPFILT(G, 'chmean', M, N, Q) Contraharmonic mean
% filtering of order Q. The
% default is Q = 1.5.
% F = SPFILT(G, 'median', M, N) Median filtering.
% F = SPFILT(G, 'max', M, N) Max filtering.
% F = SPFILT(G, 'min', M, N) Min filtering.
% F = SPFILT(G, 'midpoint', M, N) Midpoint filtering.
% F = SPFILT(G, 'atrimmed', M, N, D) Alpha-trimmed mean filtering.
% Parameter D must be a
% nonnegative even integer;
% its default value is D = 2.
%
% The default values when only G and TYPE are input are M = N = 3,
% Q = 1.5, and D = 2.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.6 $ $Date: 2003/10/27 20:07:00 $
% Process inputs.
if nargin == 2
m = 3; n = 3; Q = 1.5; d = 2;
elseif nargin == 5
Q = parameter; d = parameter;
elseif nargin == 4
Q = 1.5; d = 2;
else
error('Wrong number of inputs.');
end
% Do the filtering.
switch type
case 'amean'
w = fspecial('average', [m n]);
f = imfilter(g, w, 'replicate');
case 'gmean'
f = gmean(g, m, n);
case 'hmean'
f = harmean(g, m, n);
case 'chmean'
f = charmean(g, m, n, Q);
case 'median'
f = medfilt2(g, [m n], 'symmetric');
case 'max'
f = ordfilt2(g, m*n, ones(m, n), 'symmetric');
case 'min'
f = ordfilt2(g, 1, ones(m, n), 'symmetric');
case 'midpoint'
f1 = ordfilt2(g, 1, ones(m, n), 'symmetric');
f2 = ordfilt2(g, m*n, ones(m, n), 'symmetric');
f = imlincomb(0.5, f1, 0.5, f2);
case 'atrimmed'
if (d <= 0) | (d/2 ~= round(d/2))
error('d must be a positive, even integer.')
end
f = alphatrim(g, m, n, d);
otherwise
error('Unknown filter type.')
end
%-------------------------------------------------------------------%
function f = gmean(g, m, n)
% Implements a geometric mean filter.
inclass = class(g);
g = im2double(g);
% Disable log(0) warning.
warning off;
f = exp(imfilter(log(g), ones(m, n), 'replicate')).^(1 / m / n);
warning on;
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = harmean(g, m, n)
% Implements a harmonic mean filter.
inclass = class(g);
g = im2double(g);
f = m * n ./ imfilter(1./(g + eps),ones(m, n), 'replicate');
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = charmean(g, m, n, q)
% Implements a contraharmonic mean filter.
inclass = class(g);
g = im2double(g);
f = imfilter(g.^(q+1), ones(m, n), 'replicate');
f = f ./ (imfilter(g.^q, ones(m, n), 'replicate') + eps);
f = changeclass(inclass, f);
%-------------------------------------------------------------------%
function f = alphatrim(g, m, n, d)
% Implements an alpha-trimmed mean filter.
inclass = class(g);
g = im2double(g);
f = imfilter(g, ones(m, n), 'symmetric');
for k = 1:d/2
f = imsubtract(f, ordfilt2(g, k, ones(m, n), 'symmetric'));
end
for k = (m*n - (d/2) + 1):m*n
f = imsubtract(f, ordfilt2(g, k, ones(m, n), 'symmetric'));
end
f = f / (m*n - d);
f = changeclass(inclass, f);
⑩ 入口函數：paddedsize.m
function PQ=paddedsize(AB,CD,PARAM)
if nargin==1
PQ=2*AB;
elseif nargin==2&~ischar(CD)
PQ=AB+CD-1;
PQ=2*ceil(PQ/2);
elseif nargin==2
m=max(AB);
P=2^nextpow2(2*m);
PQ=[P,P];
elseif nargin==3
m=max([AB,CD]);
P=2^nextpow2(2*m);
PQ=[P,P];
else
error('wrong number of inputs.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [u,v]=dftuv(m,n)
u=0:(m-1);
v=0:(n-1);
idx=find(u>m/2);
u(idx)=u(idx)-m;
idy=find(v>n/2);
v(idy)=v(idy)-n;
[v,u]=meshgrid(v,u);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function g=dftfilt(f,H)
F=fft2(f,size(H,1),size(H,2));
g=real(ifft2(H.*F));
g=g(1:size(f,1),1:size(f,2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%本檔案包含paddedsize.m和dftuv.m及dftfilt.m三個檔案
f=imread('D:\test\windmill_noise.png');
figure,subplot(2,2,1),imshow(f);
F=fft2(f);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,2),imshow(s,[]); %原圖像頻譜
PQ=paddedsize(size(f));
[u,v]=dftuv(PQ(1),PQ(2));
d0=0.05*PQ(2);
F=fft2(f,PQ(1),PQ(2));
H=exp(-(u.^2+v.^2)/(2*(d0^2))); %高斯低通濾波
g=dftfilt(f,H);
subplot(2,2,3),imshow(g,[]);
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]); %高斯濾波後的圖像頻譜
d=sqrt(u.^2+v.^2);
H=1./(1+(d/d0).^2); %巴特沃茲低通濾波
g=dftfilt(f,H);
figure,subplot(2,2,1),imshow(g,[]);
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,2),imshow(s,[]); %巴特沃茲低通濾波後的圖像頻譜
h=double(d<=d0); %理想低通濾波
g=dftfilt(f,H);
H=double(d<=d0);
subplot(2,2,3),imshow(g,[]);
F=fft2(g);
fc=fftshift(F);
s=log(1+abs(fc));
subplot(2,2,4),imshow(s,[]); %理想低通濾波後的圖像頻譜