第三章频率域滤波
！！！首先本章需要定义的所有库函数代码如下 ： 一 .库函数
！！！每个代码前面有相应的题目标号注释
！！可以分题目定义也全部可以一起定义，建议本章库函数及题目代码全部放到一个文件夹中
一．库函数
%3.1定义函数
function g=gscale(f,varargin)
if length(varargin)==0
method='full8';
else method=varargin{1};
end
if strcmp(class(f),'double') & (max(f(:))>1 || min(f(:))<0)
f=mat2gray(f);
end
switch method
case 'full8'
g=im2uint8(mat2gray(double(f)));
case 'full16'
g=im2uint16(mat2gray(double(f)));
case 'minmax'
low = varargin{2};
high = varargin{3};
if low>1 || low<0 || high>1 || high<0
error('Parameters low and high must be in the range [0,1]')
end
if strcmp(class(f),'double')
low_in=min(f(:));
high_in=max(f(:));
elseif strcmp(class(f),'uint8')
low_in=double(min(f(:)))./255;
high_in=double(max(f(:)))./255;
elseif strcmp(class(f),'uint16')
low_in=double(min(f(:)))./65535;
high_in=double(max(f(:)))./65535;
end
g=imadjust(f,[low_in high_in],[low high]);
otherwise
error('Unknown method')
end
%3.1定义函数
function [out, revertclass]=tofloat(in)
%tofloat convert image to floating point
identity=@(x) x;
tosingle=@im2single;
table={'uint8', tosingle, @im2uint8
'uint16', tosingle, @im2uint16
'int16', tosingle, @im2int16
'logical', tosingle, @logical
'double', identity, identity
'single', identity, identity};
classIndex=find(strcmp(class(in),table(:,1)));
if isempty(classIndex)
error('Unsupported inut image class.')
end
out=table{classIndex,2}(in);
revertclass=table{classIndex,3};
%3.1定义函数
function [H, D] = 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 value 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.
%
% 得到指定类型的低通滤波器
% 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('Unkown filter type.')
end
%3.1定义函数
function PQ = paddedsize(AB, CD, PARAM)
% 计算填充尺寸以供基于FFT的滤波器
% PQ = PADDEDSIZE(AB),AB = [A B], PQ = 2 * AB
%
% PQ = PADDEDSIZE(AB, 'PWR2')， PQ（1） = PQ（2） = 2 ^ nextpow2(2 * m), m =
% MAX(AB).
%
% PQ = PADDEDSIZE(AB, CD)，AB = [A B], CD = [C D]
%
% PQ = PADDEDSIZE(AB, CD, 'PWR2')， PQ（1） = PQ（2） = 2 ^ nextpow2(2 * m), m =
% MAX([AB CD]).
if nargin == 1
PQ = 2 * AB;
elseif nargin == 2 & ~ischar(CD)
PQ = AB + CD -1;
PQ = 2 * ceil(PQ / 2); % ceil（N）返回比N大的最小整数，为了避免出现奇数，因为处理偶数数组快
elseif nargin == 2
m = max(AB);
P = 2 ^ nextpow2(2 * m); % nextpow2（N）返回第一个P，使得2. ^ P> = abs（N）。
% 对于FFT运算，找到最接近两个序列长度的2 的幂次方通常很有用。
PQ = [P, P];
elseif nargin == 3
m = max([AB CD]);
P = 2 ^ nextpow2(2 * m);
PQ = [P, P];
else
error('Wrong number of input')
end
%3.1 3.3定义函数
function [U, V] = dftuv(M, N)
%DFTUV Computes meshgrid frequency matrices.
% [U, V] = DFTUV(M, N) computes meshgrid frequency matrices U and
% V. U and V are useful for computing frequency-domain filter
% functions that can be used with DFTFILT. U and V are both
% M-by-N.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.3 $ $Date: 2003/04/16 22:30:34 $
% Set up range of variables.
u = 0:(M - 1);
v = 0:(N - 1);
% Compute the indices for use in meshgrid.
idx = find(u > M/2);
u(idx) = u(idx) - M;
idy = find(v > N/2);
v(idy) = v(idy) - N;
% Compute the meshgrid arrays.
[V, U] = meshgrid(v, u);
%3.2定义函数
function g = dftfilt(f, H)
% DFTFILT performs frequency domain filtering.
% G = DFTFILT(F, H) filters F in the frequency domain using the filter
% transfer function H. The output, G, is the filtered image, which has
% the same size as F. DFTFILT automatically pads F to be the same size as
% H. Function PADEDESIZE can be used to determine an appropriate size
% for H.
% G = DFTFILT（F，H）使用滤波器传递函数H在频域中对输入图像F滤波。 输出G是滤波后的图像，
% 其大小与F相同。DFTFILT自动将F填充到与H相同的大小 ，PADEDESIZE函数可用于确定H的合适大小。
%
% DFTFILT assumes that F is real and that H is a real, uncentered,
% circularly-symmetric filter function.
% DFTFILT假设F是实数，H是实数，未中心，循环对称的滤波函数。
% Obtain the FFT of the padded input.
% 获取填充之后的FFT变换
F = fft2(f, size(H, 1), size(H, 2));
% Perform filtering
% 滤波
g = real(ifft2(H .* F));
% Crop to orihinal size
% 剪切到原始尺寸
g = g(1 : size(f, 1), 1 : size(f, 2));
%3.5 3.6 定义函数
function H = hpfilter(type, M, N, D0, n)
% LPFILTER 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). To view the
% filter as an image or mesh plot, it should be centered using H =
% fftshift(H)
% Valid value 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.
% The transfer function Hhp of 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.
%
% 计算给定类型（理想、巴特沃兹、高斯）的频域高通滤波器
% Use function dftuv to set up the meshgrid arrays needed for computing the
% required distances.
if nargin == 4
n = 1;
end
Hlp = lpfilter(type, M, N, D0, n);
H = 1 - Hlp;
二.课本例题
3.1 习题代码及结果图
clc
clear
f = imread('F:\shi yan\图像处理\1.tif');
[M,N]=size(f);
[f,revertclass]=tofloat(f);
F=fft2(f);
sig=10;
H=lpfilter('gaussian',M,N,sig);
G=H.*F;
g=ifft2(G);
g=revertclass(g);
PQ = paddedsize(size(f)); %首先进行尺寸填充
Fp = fft2(f,PQ(1), PQ(2)); %以填充后的尺寸进行FFT
Hp = lpfilter('gaussian', PQ(1), PQ(2), 2*sig); %因为现在滤波器尺寸为未填充时候的2倍，所以选用2*sig为参数
Gp = Hp.*Fp; %让频谱和滤波器相乘得到Gp
gp = real(ifft2(Gp)); %对Gp取IFFT后取实部
gpc = gp(1:size(f,1), 1:size(f,2)); %裁剪尺寸为原始大小
h =fspecial('gaussian',15,7);
gs = imfilter(f,h);
subplot(2,2,1);imshow(f);title('原图像');
subplot(2,2,2);imshow(g);title('巴特沃兹高通滤波');
subplot(2,2,3);imshow(gpc);title('使用填充的频域低通滤波处理后的图像') ;
subplot(2,2,4);imshow(gp);title('使用填充的频域低通滤波处理后的图像') ;
PQ = paddedsize(size(f));%使用paddedsize获得填充参数；如果输入是彩色图像，必须要灰度化rgb2gray
Hp = lpfilter('gaussian', PQ(1), PQ(2), 2*sig);
Fp = fft2(f, PQ(1), PQ(2));%得到使用填充的傅里叶变化
%生成一个大小为PQ(1) X PQ(2) 的滤波函数H。如果该滤波函数已居中，使用前要令H = fftshift(H)
Gp = Hp.*Fp;%将变换乘以滤波函数
gp = real(ifft2(G));%获得G的傅里叶逆变换的实部
gpc = gp(1:size(f,1), 1:size(f, 2));%将左上部分的矩形剪切为原来尺寸大小
gpc=revertclass(gpc);
h=fspecial('gaussian',15,7);
gs=imfilter(f,h);
实验结果图：
3.2习题代码及结果图
%3.2
f = imread('F:\shi yan\图像处理\2.tif');
%获得图像f的傅里叶频谱
F = fft2(f);
S = fftshift(log(1 + abs(F)));
S = gscale(S);
%生成空间滤波器
h = fspecial('sobel')';
%查看相应频域滤波器的图形
freqz2(h)
%滤波器本身通过如下命令获得
PQ = paddedsize(size(f));
H = freqz2(h, PQ(1), PQ(2));
H1 = ifftshift(H);
%显示H和H1的绝对值
%在空间域中滤波后图像
%该语句默认用0填充边界。
gs = imfilter(double(f), h);
%频域滤波后图像
gf = dftfilt(f, H1);
d = abs(gs - gf);
max(d(:))
min(d(:))%通过计算最大差和最小差说明
%空间滤波和频域滤波得到的图像相同
subplot(4, 3, 1), imshow(f);
subplot(4, 3, 2), imshow(S);
subplot(4, 3, 3), imshow(abs(H), []);
subplot(4, 3, 4), imshow(abs(H1), []);
subplot(4, 3, 5), imshow(gs, []);
subplot(4, 3, 6), imshow(gf, []);
subplot(4, 3, 7), imshow(abs(gs), []);
subplot(4, 3, 8), imshow(abs(gf), []);
subplot(4, 3, 9), imshow(abs(gs > 0.2*abs(max(gs(:)))));
subplot(4, 3, 10), imshow(abs(gf > 0.2*abs(max(gf(:)))));
实验结果图：
3.3习题代码及结果图
%三步走输入
[U,V]=dftuv(8,5);
DSQ=U.^2+V.^2;
fftshift(DSQ)
实验结果图：
3.4习题代码及结果图
%3.4
%高斯低通
I=imread('F:\shi yan\图像处理\4.tif');
subplot(221),imshow(I);
title('原图像');
Y=fft2(im2double(I));%傅里叶变换
Y=fftshift(Y);%频谱搬移,直流分量搬移到频谱中心
subplot(222), imshow(log(abs(Y)+1),[]);
title('图像傅里叶变换取对数所得频谱');
[M,N]=size(Y);%获得图像的高度和宽度
h=zeros(M,N);%滤波器函数
%图像中心点
M0=M/2;
N0=N/2;
%截至频率距离圆点的距离，delta表示高斯曲线的扩散程度
D0=40;
delta=D0;
for x=1:M
for y=1:N
%计算点（x,y）到中心点的距离
d2=(x-M0)^2+(y-N0)^2;
%计算高斯滤波器
h(x,y)=exp(-d2/(2*delta^2));
end
end
%滤波后结果
res=h.*Y;
res=real(ifft2(ifftshift(res)));
subplot(223),imshow(res);
title('高斯低通滤波所得图像');
subplot(224),imshow(h);
title('高斯低通滤波器图象');
实验结果图：
3.5习题代码及结果图
%3.5
H = fftshift(hpfilter('gaussian', 500, 500, 50));
mesh(double(H(1:10:500,1:10:500)));
axis tight;
subplot(3, 2, 1), surf(double(H(1:10:500, 1:10:500))), title('1');
B = fftshift(hpfilter('gaussian', 500, 500, 50));
mesh(double(B(1:10:500,1:10:500)));
axis tight;
colormap([0 0 0]);
axis off;
subplot(3, 2, 2), surf(double(B(1:10:500, 1:10:500))), title('2');
G = fftshift(hpfilter('gaussian', 500, 500, 50));
mesh(double(G(1:10:500,1:10:500)));
axis tight;
colormap([0 0 0]);
axis off;
view(-25,30);
subplot(3, 2, 3), surf(double(G(1:10:500, 1:10:500))), title('3');
K = fftshift(hpfilter('gaussian', 500, 500, 50));
mesh(double(K(1:10:500,1:10:500)));
axis tight;
colormap([0 0 0]);
axis off;
view(-25,0);
subplot(3, 2,4), surf(double(K(1:10:500, 1:10:500))), title('4');
实验结果图：
3.6习题代码及结果图
%3.6
H = fftshift(hpfilter('ideal', 500, 500, 50));
mesh(H(1:10:500,1:10:500));
axis([0 50 0 50 0 1]);
colormap([0 0 0])
axis off
grid off
B = fftshift(hpfilter('btw', 500, 500, 50));
mesh(B(1:10:500,1:10:500));
axis([0 50 0 50 0 1]);
colormap([0 0 0]);
axis off;
grid off;
G = fftshift(hpfilter('gaussian', 500, 500, 50));
mesh(G(1:10:500,1:10:500));
axis([0 50 0 50 0 1]);
colormap([0 0 0]);
axis off;
grid off;
subplot(3, 2, 1), surf(double(H(1:10:500, 1:10:500))), title('理想高通滤波器绘图');
subplot(3, 2, 2), imshow(H,[]), title('理想高通滤波器图像');
subplot(3, 2, 3), surf(double(B(1:10:500, 1:10:500))), title('布特沃斯高通滤波器透视图');
subplot(3, 2, 4), imshow(B,[]), title('布特沃斯高通滤波器图像');
subplot(3, 2, 5), surf(double(G(1:10:500, 1:10:500))), title('高斯高通滤波器透视图');
subplot(3, 2, 6), imshow(G,[]), title('高斯高通滤波器图像');
实验结果图：
3.7习题代码及结果图
%3.7
clear
d0=50; %阈值
image=imread('F:\shi yan\图像处理\7.tif');
[M ,N]=size(image);
img_f = fft2(double(image));%傅里叶变换得到频谱
img_f=fftshift(img_f); %移到中间
m_mid=floor(M/2);%中心点坐标
n_mid=floor(N/2);
h = zeros(M,N);%高斯低通滤波器构造
for i = 1:M
for j = 1:N
d = ((i-m_mid)^2+(j-n_mid)^2);
h(i,j) = 1-exp(-(d)/(2*(d0^2)));
end
end
img_lpf = h.*img_f;
img_lpf=ifftshift(img_lpf); %中心平移回原来状态
img_lpf=uint8(real(ifft2(img_lpf))); %反傅里叶变换,取实数部分
subplot(1,2,1);imshow(image);title('原图');
subplot(1,2,2);imshow(img_lpf);title('高斯高通滤波');
实验结果图：
3.8习题代码及实验结果图
%3.8
clc
clear
f=imread('F:\shi yan\图像处理\8.tif'); %读入图像
subplot(2,2,1);
imshow(f)
title('原始图像')
%高斯高通滤波
I=double(f);
g=fft2(I);
g=fftshift(g);
[M,N]=size(g);
D0=40;
m=fix(M/2);n=fix(N/2);
for i=1:M
for j=1:N
D=sqrt((i-m)^2+(j-n)^2);
H=exp(-(D.^2)./(2*(D0^2)));
result(i,j)=(1-H)*g(i,j);
end
end
result=ifftshift(result);
J1=ifft2(result);
J2=uint8(real(J1));
subplot(2,2,2);
imshow(J2)
title('高斯高通滤波后的图像')
%高频强调滤波
F=0.5+0.75*(1-H);
G=F*g;
result2=ifftshift(G);
J3=ifft2(result2);
J4=uint8(real(J3));
subplot(2,2,3)
imshow(J4)
title('高频强调滤波后的图像')
%对高频强调滤波后图像进行直方图均衡化
J5=histeq(J4,256);
J5=uint8(J5);
subplot(2,2,4);
imshow(J5)
title('直方图均衡化后的图像')
实验结果图：
这样第三章的课本习题就完成了
三、小编提示！
实验相关图片可去数字图像处理课本官网下载压缩包！网站：https://imageprocessingplace.com/DIP-3E/dip3e_book_images_downloads.htm，保存的实验图片记得进行修改名称以便代码读取图片喔！
相关阅读
发表评论