三、基于MATLAB的数字图像处理————频率域滤波(一)

(三)基于MATLAB的数字图像处理————频率域滤波

一、二维离散傅里叶变换
令f(x,y)表示一幅大小为M*N像素的图像,其中x=0,1,2,…,M-1,y=0,1,2,…N-1
由**F(u,v)*表示的二维离散傅里叶变换(DFT)为:
在这里插入图片描述
式中,u=0,1,2…,M-1,v=0,1,2,…,N-1。显然,F也是大小为MN的矩形区域,通常称为频率矩形。
离散傅里叶反变换(IDFT)的形式为:
在这里插入图片描述
式中,x=0,1,2,…,M-1
y=0,1,2,…,N-1
。由 F(u,v) 可以得到 f(x,y)
对于 f(x,y) 为实函数,他的 DFT 变换具有周期性:
在这里插入图片描述
他的 IDFT 也具有周期性:
在这里插入图片描述
R(u,v)I(u,v) 分别表示 F(u,v) 的实部和虚部,则傅里叶谱定义为
在这里插入图片描述
变换的相角定义为
在这里插入图片描述
F(u,v) 还能表示为:
在这里插入图片描述
二、观察DFT的相关函数
在MATLAB中,一个图像数组的DFT可以用函数 fft2 得到

F=fft2(f)

若需要进行零填充,得到大小为P*Q的结果则有语法:

F=fft2(f,P,Q)
S=abs(F)           %得到F的傅里叶谱
Fc=fftshift(F)     %将变换的原点移动到频率矩形的中心
F=ifftshift(Fc)    %将居中的结果反转
phi=angle(F)       %得到DFT的相角
phi=atan2(imag(F),real(F))      %同样得到DFT的相角
F=S.*exp(i*phi)    %由 给出的傅里叶谱和相角得到DFT

因为 fft2 把输入图像毫无缩放地变成了double类,但它值的范围仍然是 [0 255] 这样在进行图像显示时就会出错,我们定义 tofloat 函数将输入图像转化为范围为 [0 1] 内的浮点数,然后在过程的末尾用 tofloatrevertclass 特性将结果转换成与初始图像相同的类。代码如下:

function [out,revertclass]=tofloat(in)
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 input image class.');
end
out=table{classindex,2}(in);
revertclass=table{classindex,3};

三、频率域滤波
1,基本原理
空间域滤波和频率域滤波的基础都是卷积,卷积定理表述为:空间域乘积的DFT为频域卷积;空间域卷积的DFT为频域内的乘积。
在进行滤波处理时我们经常用到的是后者:
在这里插入图片描述
H(u,v) 为滤波器传递函数,频域滤波的思想就是选择一个滤波器传递函数去对原图像的DFT进行处理(相乘),将结果再IDFT,得到的即为滤波之后的结果。
为了进行DFT运算,构造一个将大小为AD的f(x,y)和CD的h(x,y)构造成大小均为P*Q的扩展函数,要求
在这里插入图片描述

在这里插入图片描述
实现的函数代码:

function PQ=paddedsize(AB,CD,PARAM)
if nargin == 1
    PQ=2*AB;
elseif nargin == 2 &&~ischar(CD) %CD为字符串返回真
    PQ=AB+CO-1;
    PQ=2*ceil(PQ/2);  %ceil(x)返回大于等于x的最小整数
elseif nargin==2
    m=max(AB);
    p=2^nextpow2(2*m); %nextpow2(n)返回大于或等于n的绝对值的2的最小整数次幂
    PQ=[P,P];
elseif (nargin==3) && strcmpi(PARAM,'pwr2')
    m=max([AB CD]);
    P=2^nextpow2(2*m);
    PQ=[P,P];
else
    error('wrong number of inputs.')
end

2,基本步骤
a,使用tofloat函数把输入图像转换为浮点图像

[f,revertclass]=tolfloat(f);

b,使用函数paddedzsize获得填充参数

PQ=paddedzsize;

c,得到的有填充图像的傅里叶变换

F=fft2(f,PQ(1),PQ(2));

d,生成一个大小为PQ(1)*PQ(2)的滤波器函数。(如果滤波器函数是居中的,则在使用滤波器之前要令H=ifftshift(H))。

e,用滤波器乘以该变换

G=H.*F;

f,获得G的IFFT:

g=ifft2(G);

g,将左上部的矩形修剪为原始大小:

g=g(1:size(f,1),1:size(f,2));

h,需要时,将滤波后的图像转换成原始输入图像的类:

g=revertclass(g);

3,用于频域滤波的M函数
假定滤波器函数的大小已被适当的调整,这时,我们将滤波后的图像转换成与输入想同的类,或者我们需要处理浮点型结果。定义一个M函数:

function g=dftfilt(f,H,classout)

[f,revertclass]=tofloat(f);
F=fft2(f,size(H,1),size(H,2));
g=ifft2(H.*F);
g=g(1:size(f,1),1:size(f,2));
if nargin==2||strcmp(classout,'original')
    g=revertclass(g);
elseif strcmp(classout,'fltpoint')
    return
else
    error('undefined class for the output image')
end

4,从空间滤波器获得频率域滤波器
图像处理函数 freqz2可以实现从空间滤波器获得频率域滤波器的转换。

H=freqz2(h,R,C)

其中,h是一个二维空间滤波器,H是相应的二维频率滤波器。R是行数,C 是我们希望的滤波器的列数。通常,我们令R=PQ(1),C=PQ(2)。

做一个空间滤波和频率域滤波的比较

f=tofloat(f);
>> F=fft2(f);
>> S=fftshift(log(1+abs(F)));
>> imshow(S,[])
>> h=fspecial('sobel')
>h =

     1     2     1
     0     0     0
    -1    -2    -1
>> PQ=paddedsize(size(f));
>> H=freqz2(h,PQ(1),PQ(2));
>> H1=ifftshift(H);
>> imshow(abs(H),[])
>> figure,imshow(abs(H1),[])

原始图像以及傅里叶谱为
灰度图像及其傅里叶谱
由空间滤波器得到的频率域滤波器的傅里叶谱及其ifftshift处理后的傅里叶谱
在这里插入图片描述

进行空间滤波

gs=imfilter(f,h);
imshow(gs,[])

在这里插入图片描述
进行频率域滤波

gf=dftfilt(f,H1);
 imshow(gf,[])

在这里插入图片描述
图像中的灰色调是由于gs和gf的负数值引起的,通过标定命令imshow,负数值会使得平均值增大。上面生成的Sobel模板h通过使用响应的绝对值来检测图像的水平边缘,更有意义。

imshow(abs(gs),[])
 figure,imshow(abs(gf),[])

空间滤波和频率域滤波的绝对值结果对比
在这里插入图片描述

使用空间滤波和频率域滤波得到的图像所有实用目的是相同的,通过计算他们的差值,我们证明这一事实:

 d=abs(gs-gf);
>> max(d(:))

ans =

  single

  8.0466e-07  %在这里 这个值的大小可以忽略不计。

>> min(d(:))

ans =

  single

     0

这个例子说明对于一幅图像的滤波,无论是选择空间滤波,还是选择这种空间滤波相应的频域滤波,最后得到的结果都是相同的。

  • 16
    点赞
  • 88
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值