数字图像处理第三章

最新推荐文章于 2022-10-02 11:22:53 发布

Daos2020

最新推荐文章于 2022-10-02 11:22:53 发布

阅读量2.4k

点赞数

文章标签：计算机视觉

本文链接：https://blog.csdn.net/weixin_44086996/article/details/108591774

版权

二维离散傅里叶变换

（一）理论推导
（二）在 MATLAB 中计算及观察二维 DFT
（三）频域滤波

背景知识：傅里叶变换在时域和频域上都呈离散的形式，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作经过周期延拓成为周期信号再作变换。

（一）理论推导

令f(x,y)代表一幅大小为 MxN的数字图像，其中x=0,1,2,…,M-1, y=0,1,2,…,N-1, 由 F(u,v)表示的f(x,y)的二维离散傅立变换(DFT)可以由下式给出:
在这里插入图片描述
其中u=0,1,2,…,M-1,v=0,1,2,…,N-1。我们可以借助决定频率的变量 u、v(x和 y求和) 将指数形式展开为正弦和余弦的形式。频域系统是以 u,v(频率)为变量来表示 F(u, v)的坐标系。频域矩形和输入图像是同等大小的。离散傅立叶反变换(IDFT)的形式为：
在这里插入图片描述
其中x=0,1,2,…,M-1,y=0,1,2,…,N-1。因此，给定 F(u， v), 我们可以借助 IDFT(傅立叶反变换)来得到f(x,y)。

特别地，由于 MATLAB 中的数组索引以 1 开头，而不是以 0 开头，MATLAB 中的F(1,1)和f(1,1)对应数学公式中正变换和反变换的 F(0,0)和f(0,0)。

傅立叶谱可以定义如下：
在这里插入图片描述
注：傅里叶频谱是一种理想的可用于描绘周期或者近似周期的二维图像模式的方向性的方法。而频谱法正是基于傅里叶频谱的一种纹理描述方法。全局纹理模式在空域中很难检测出来，但是转换到频域中则很容易分辨。因此，频谱纹理对区分周期模式或非周期模式以及周期模式之间的不同十分有效。通常，全局纹理模式对应于傅里叶频谱中能量十分集中的区域，即峰值突起处。

变换的相角定义如下:
在这里插入图片描述
这两个函数可在极坐标形式下用于表示复函数F(u,v)：

功率谱可以定义为幅度的平方:

在这里插入图片描述
如果f(x,y)是实函数，它的傅立叶变换就是关于原点共轭对称的:

在这里插入图片描述
这意味着傅立叶谱也是关于原点对称的:

在这里插入图片描述
可以将它直接代入到公式 F(u,v)中:

在这里插入图片描述
其中， $k_1$ 和 $k_2$ 是整数。DFT 在u, v 方向上是无穷周期的，周期由 M和N决定。周期性也是 DFT 逆变换的重要特性：

在这里插入图片描述
也就是说，通过傅立叶反变换得到的图像也是无穷周期的。并且，在计算离散傅立叶变换时只计算它的一个周期，即仅处理尺寸为MXN的数组。

（二）在 MATLAB 中计算及观察二维 DFT

在实践中，DFT 及其反变换可以用快速傅立叶变换(FFT)算法实现。一幅图像数组 f 的 FFT 可以在 MATLAB中用函数 fft2 得到，下面先引入一段示例代码及其运行结果，在对其分析。
示例代码：

>> I = rgb2gray(imread('test914.png'));
>> F1 = fft2(I);
>> F2 = fft2(I,200,200);
>> F2c = fftshift(F2);
>> S2 = log(1+abs(F2c));
>>> Fi = ifftshift(abs(F2c));
>> subplot(231),imshow(I,[]);title('a')
>> subplot(232),imshow(F1,[]);title('b')
>> subplot(233),imshow(F2,[]);title('c')
>> subplot(234),imshow(F2c,[]);title('d')
>> subplot(235),imshow(S2,[]);title('e')
>> subplot(236),imshow(Fi,[]);title('f')
>>

运行结果：
在这里插入图片描述
对实验分析：

图a到b是将原图进行傅里叶变换得到的频谱，可以看到图b的四周的较亮区域，说明该处的频率大，变化剧烈。
为了便于观察，将远点移至中心(fftshift)。移动之后如图d所示，可以清晰的看到图像中心的频率交大。
图b到c的操作是0填充的操作，这是因为使用傅里叶变换滤波时，需要对数据进行0填充，此处暂时没有用到。
为了观测的更清晰，中心处的亮度值占支配地位，可以通过log变换解决这个问题。图e显示的是使用log函数均衡化的结果。

我们使用另外一张图像观察下傅里叶频谱

I = rgb2gray(imread('heibai.jpg'));
>> F2 = fft2(I,200,200);
>> F2c = fftshift(F2);
>> S2 = log(1+abs(F2c));
>>> figure;subplot(121),imshow(I,[]);title('a')
>> subplot(122),imshow(S2,[]);title('b')

在这里插入图片描述
从图中我们可以看出傅里叶频谱高亮部分正是原图变化剧烈的地方。实验的结果正如预期的一样。下面进入第三部分频域滤波

（三）频域滤波

3.1 基础知识

空(间)域和频(率)域的线性滤波的基础都是卷积定理，可以被写作:
在这里插入图片描述
逆变换为:

在这里，符号"星号"表示两个函数的卷积，双箭头两边的表达式组成傅立叶变换对。例如, 第一个表达式表明两个空间函数(表达式左侧的项)的卷积可以通过计算两个傅立叶变换函数(表达式的右侧)的乘积的反变换得到。相反，两个空间函数的卷积的正傅立叶变换给出了两个函数傅立叶变换的乘积。在处理频域滤波的时候，输入的傅立叶变换不必居中处理。

当工作在离散量时，我们知道F和H是周期性的，这意味着在离散频域执行卷积也是周期性的。由于这个原因，用DFT 执行的卷积叫做循环卷积。为了保证空间和循环卷积给出相同结果的，所以在滤波操作前要先填充适当的0。注：当乘以居中处理后的函数 F(u，v)后，衰减F(u，v)的高频分量,相应的低频分量没有变化。具有这种特性的滤波器叫做低通滤波器。

补0的方法如下：

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 input.')

end

上述语法对尺寸是 PQ⑴xPQ(2)的结果图像添加了足够多的 0。这是为了避免折叠误差。

下面我们使用高斯低通滤波器处理图像。首先我们观察一下波形器的图像：

lpfilter('gaussian',M,N,sig);

在这里插入图片描述
其中滤波器的大小与待处理图像相同（大约800*1300），sig = 10。可以看出该滤波器对高频的分量处理，相对的低频部分没有什么变化。下面使用滤波器处理图像的示例代码如下：

>> f= rgb2gray(imread('dog.png'));
>> clear
>> f= rgb2gray(imread('dog.png'));
>> [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));
>> HP = lpfilter('gaussian',PQ(1),PQ(2),2*sig);
>> GP = HP.*FP;
>> gp = ifft2(GP);
>> gpc = gp(1:size(f,1),1:size(f,2));
>> gpc = revertclass(gpc);
>>> subplot(131),imshow(f),title('原图');
>> subplot(132),imshow(g),title('未填充');
>> subplot(133),imshow(gpc,[]),title('填充');

在这里插入图片描述
可以看出填充与未填充都将原图像模糊了，但可以观察到填充图像与未填充图像的边缘差别较大，这是因为填充之后当一块亮区和一块暗区处在滤波器之下时，图像会变灰，模糊了输出。

3.2DFT 滤波的基本步骤

前面的讨论可以概括为下面几个步骤，其中 f 是被滤波处理的图像，gpc为处理结果，同时假设滤波函数 H与填充后的图像大小相等。

1 把输入图像变成double类型图像：

f = im2double(rgb2gray(h));

2用函数 paddedsize 获得填充参数：

PQ = paddedsize(size(f));

3 得到有填充的傅立叶变换:

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

4用滤波器乘以FFT 变换：

Gp=Hp.*Fp;

5 获得 Gp 的逆 FFT 变换:

gp=ifft2(Gp);

6 修剪左上部矩形为原始大小:

gpc = gp(1:size(f,1), 1:size(f,2));

7 把滤波过的图像变换为输入图像的类，如果需要。

g = revertclass(g);

3.3从空域滤波器获得频域滤波器

背景知识：空间域、频率域、时域。
1，空间域：

空间域（spatial domain）也叫空域，即所说的像素域，在空域的处理就是在像素级的处理，如在像素级的图像叠加。通过傅立叶变换后，得到的是图像的频谱。表示图像的能量梯度。

2，频率域：

   频率域（frequency domain。）任何一个波形都可以分解成多个正弦波之和。每个正弦波都有自己的频率和振幅。所以任意一个波形信号有自己的频率和振幅的集合。频率域就是空间域经过傅立叶变换的信号

3，时域：

   时域（时间域）——自变量是时间,即横轴是时间,纵轴是信号的变化。其动态信号x（t）是描述信号在不同时刻取值的函数。

傅里叶变换是实现从空域或时域到频域的转换工具。

空域滤波器

使用空域模板进行的图像处理，被称为空域滤波。模板本身被称为空域滤波器。

首先使用空域滤波器对图像处理，生成一个空域滤波器（‘sobel’）。示例代码：

>> h = fspecial('sobel')'

h =

     1     0    -1
     2     0    -2
     1     0    -1

我们可以观察相应频域滤波器的图形

>>freqz2(h)

在这里插入图片描述
接下来设置一个频域滤波器：

>> PQ = paddedsize(size(f));
>> H = freqz2(h,PQ(1),PQ(2));
>> H1 = ifftshift(H);

这里扩充是必须的，使得原点在频域矩形的左上角。
分别使用空域与频域滤波器对图像处理：

gs = imfilter(f,h);     %空域
gf = dftfilt(f,H1);
>> figure, imshow(abs(gf) > 0.2*abs(max(gf(:))))
>> figure, imshow(abs(gs) > 0•2*abs(max(gs(:))))
>>> figure;imshow(f)

在这里插入图片描述

代码运行效果显示空域处理与频域处理图像显示一样，我们再通过计算它们的差别来证实这个事实：

d = abs(gs - gf);

>> max(d{:))
ans
1.2973e-006

结论：可以看出用空域和频域滤波得到的图像对于实际应用的目的是一样的。

注：dftfilt为频域滤波的 M-函数，其方法如下：

function g = dftfilt(f, H)
%   G = DFTFILT（F，H）使用滤波器传递函数H在频域中对输入图像F滤波。 输出G是滤波后的图像，
%   其大小与F相同。DFTFILT自动将F填充到与H相同的大小 ，PADEDESIZE函数可用于确定H的合适大小。
%   DFTFILT假设F是实数，H是实数，未中心，循环对称的滤波函数。
% 获取填充之后的FFT变换
F = fft2(f, size(H, 1), size(H, 2));
g = real(ifft2(H .* F));
% 剪切到原始尺寸
g = g(1 : size(f, 1), 1 : size(f, 2));

3.4在频域中直接生成滤波器

3.4.1 建立网格数组以实现频域滤波器

对于 M-函数最主要的是：需要在频率矩形中计算任意点到规定点的距离函数。称为dftuv的M-函数提供了距离计算以及其他类似应用所需要的网格数组。由 dftuv 生成的网格数组已经满足 fft2 和 ifft2 处理的需要，不需要数据重排。
dftuv方法如下：

function [U, V] = dftuv(M, N)
% [U，V] = DFTUV（M，N）计算网格频率矩阵U和V。 U和V对于计算可与DFTFILT一起使用的
% 频域滤波器函数很有用。 U和V都是M-by-N。

% 设置变量范围
u = 0 : (M - 1);
v = 0 : (N - 1);
% 计算网格的索引，即将网络的原点转移到左上角，因为FFT计算时变换的原点在左上角。
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);

使用dftuv生成网格，运行结果如下：

>> [U,V] = dftuv(8,5);
>> DSQ = U.^2 + V.^2

DSQ =

     0     1     4     4     1
     1     2     5     5     2
     4     5     8     8     5
     9    10    13    13    10
    16    17    20    20    17
     9    10    13    13    10
     4     5     8     8     5
     1     2     5     5     2

此时距离在左上角是 0,最大距离位置在频域矩形中心，可以使用函数 fftshift 来获得关于频域矩形中心的距离：

>> fftshift(DSQ)

ans =

    20    17    16    17    20
    13    10     9    10    13
     8     5     4     5     8
     5     2     1     2     5
     4     1     0     1     4
     5     2     1     2     5
     8     5     4     5     8
    13    10     9    10    13

可以看到原点坐标移动到了中心。函数 hypot 执行与 D = sqrt ( U.^ 2 + V.^2) 相同的计算.

>> [U,V] = dftuv(8,5);
>> DSQ = hypot(U,V)

DSQ =

         0    1.0000    2.0000    2.0000    1.0000
    1.0000    1.4142    2.2361    2.2361    1.4142
    2.0000    2.2361    2.8284    2.8284    2.2361
    3.0000    3.1623    3.6056    3.6056    3.1623
    4.0000    4.1231    4.4721    4.4721    4.1231
    3.0000    3.1623    3.6056    3.6056    3.1623
    2.0000    2.2361    2.8284    2.8284    2.2361
    1.0000    1.4142    2.2361    2.2361    1.4142

>> fftshift(DSQ)

ans =

    4.4721    4.1231    4.0000    4.1231    4.4721
    3.6056    3.1623    3.0000    3.1623    3.6056
    2.8284    2.2361    2.0000    2.2361    2.8284
    2.2361    1.4142    1.0000    1.4142    2.2361
    2.0000    1.0000         0    1.0000    2.0000
    2.2361    1.4142    1.0000    1.4142    2.2361
    2.8284    2.2361    2.0000    2.2361    2.8284
    3.6056    3.1623    3.0000    3.1623    3.6056

>>

3.4.2频域低通(平滑)滤波器
理想低通滤波器(ILPF)具有如下传递函数：
在这里插入图片描述
其中， $D_0$ 为正数，D(u,v)为点(u,v)到滤波器中心的距离。满足的点的轨迹为圆。因为滤波器H(u,v)乘以一幅图像的傅立叶变换，我们看到理想滤波器切断(乘以0)圆外的所有 F(u,v)分量，而保留圆上和圆内的点不变(乘以 1)。虽然这个滤波器不能用类似的电子元件实现, 但的确可以在计算机中用前面介绍的传递函数实现。理想滤波器的特性通常用来解释振铃和折叠误差等现象。图像处理中，对一幅图像进行滤波处理，若选用的频域滤波器具有陡峭的变化，则会使滤波图像产生“振铃”，所谓“振铃”，就是指输出图像的灰度剧烈变化处产生的震荡，就好像钟被敲击后产生的空气震荡。

下面我们根据上式子编写理想低通滤波器：

>> f = (imread('hm916.jpg'));
>> g = fftshift(fft2(f));
>> [U,V] = size(f);
>> D0 = 20;      
>> n1 = fix(U/2);
>> n2 = fix(V/2);     % n1，n2为中心坐标
>> for i = 1:U;
for j = 1:V;
d = sqrt((i-n1)^2+(j-n2)^2);   %计算各点距离中心的点
if d>D0
h = 0;
else
h = 1;
end
result(i,j) = h*g(i,j);
end
end
>> result = ifftshift(result);
>> M = ifft2(result);
>> subplot(121),imshow('hm916.jpg'),title('原图');
>> subplot(122),imshow(real(M)),title('滤波之后');

在这里插入图片描述
n阶的布特沃斯低通滤波器(BLPF), 具有从滤波器中心到 A) 的距离的截止频率，传递函数为:

与理想低通滤波器不同，布特沃斯低通滤波器的传递函数在 A) 点没有尖锐的不连续。

现在对图像进行布特沃斯低通滤波：

>> I1=imread('b.jpg');
f=double(I1);                %数据类型转换
g=fft2(f);                   %图像傅里叶转换
g=fftshift(g);               %傅里叶变换平移
F2=log(abs(g));              %对傅里叶变换结果取绝对值，然后取对数
[N1,N2]=size(g);             %傅里叶变换图像尺寸
n=2;                         %参数赋初始值
d0=5;
n1=fix(N1/2);                %数据圆整
n2=fix(N2/2);                %数据圆整
for i=1:N1                   %遍历图像像素
for j=1:N2 
d=sqrt((i-n1)^2+(j-n2)^2);
if d==0
h=0;
else
h=1/(1+(d/d0)^(2*n));
end
result(i,j)=h*g(i,j);         %图像矩阵计算处理
end
end
result=ifftshift(result);
X2=ifft2(result);
X3=uint8(real(X2));
subplot(1,2,1),imshow(I1),title('原图像');
subplot(1,2,2),imshow(X3),title('Butterworth低通滤波图像');

在这里插入图片描述
对比图像使用理想低通滤波器与布特沃斯低通滤波器，可以看出布特沃斯滤波之后的图像更顺滑，是因为布特沃斯低通滤波器的传递函数在 A) 点没有尖锐的不连续。

滤波器就是建立的一个数学模型，通过这个模型来将图像数据进行能量转化，噪声就是属于高频率部分，高斯滤波器平滑处理后降低噪声的影响。若使用理想滤波器，会在图像中产生振铃现象。采用高斯滤波器的话，系统函数是平滑的，避免了振铃现象。它的特性是连续性衰减，而不像理想滤波器那样陡峭变化，即明显的不连续性。因此采用该滤波器滤波在抑制噪声的同时，图像边缘的模糊程度大大减小，没有振铃效应产生。可用于平滑处理，如图像由于量化不足产生虚假轮廓时，常可用低通滤波进行平滑以改进图像质量。传递函数由下式给出：
在这里插入图片描述

现在对图像进行高斯低通滤波：

>> w = imread('g.jpg');
f = im2double(rgb2gray(w));    %把图像变为灰度图像
PQ = paddedsize(size(f));      %产生滤波时所需大小的矩阵
[U, V] = dftuv(PQ(1),PQ(2)); 
DO = 0.05*PQ(2);               %设定高斯滤波器的阈值
F = fft2(f, PQ(1), PQ(2));     %傅里叶变换
H = exp(-(U.^2+V.^2)/(2*(DO^2))); 
g = dftfilt(f,H);              %用滤波器对图像进行频域滤波
figure;
subplot(2, 2, 1), imshow(w), title('原图像');
subplot(2, 2, 2), imshow(fftshift(H),[]), title('以图像显示的高斯低通滤波器');
subplot(2, 2, 3), imshow(log(1 + abs(fftshift(F))),[]), title('对数增强并居中的fft');
subplot(2, 2, 4), imshow(g,[]), title('高斯低通滤波后的图像的谱');

在这里插入图片描述

3.4.3线框及表面绘制
使用函数 mesh绘制线框图：

H = fftshift(lpfilter('gaussian', 500, 500, 50)); 
mesh(double(H(1:10:500, 1:10:500)));
axis tight;

运行效果如下：
在这里插入图片描述
线框图默认为彩色的，从底部为蓝色渐变到顶部为红色。通过键入下面的命令，可以把绘图的线条变为黑色，并消除轴和栅格：

H = fftshift(lpfilter('gaussian', 500, 500, 50)); 
mesh(double(H(1:10:500, 1:10:500)));
colormap ( [0 0 0] );   %把绘图的线条变为黑色，
axis off ;      %消除轴
grid off ;      %消除栅格

在这里插入图片描述
以表面图代替线框图的绘制函数：

H = fftshift(lpfilter('gaussian', 500, 500, 50)); 
surf(double(H(1:10:500, 1:10:500)));
axis tight;
colormap(gray);   %将颜色转为灰度
axis off;

在这里插入图片描述
使用mesh函数绘制上一小节的理想低通滤波器、布特沃斯低通滤波器、高斯低通滤波器：

>> d0=8;
M=60;N=60;
c1=floor(M/2);     
c2=floor(N/2);      
h1=zeros(M,N);      %理想低通滤波
h2=zeros(M,N);      %布特沃斯滤波
h3=zeros(M,N);      %高斯滤波
sigma=4;
n=4;%布特沃斯阶数
for i=1:M
    for j=1:N
        d=sqrt((i-c1)^2+(j-c2)^2);
        if d<=d0
            h1(i,j)=1;
        else
            h1(i,j)=0;
        end
        h2(i,j)=1/(1+(d/d0)^(2*n)); 
        h3(i,j)=exp(-d^2/(2*sigma^2)); 
    end
end
 mesh(double(h1(1:10:60, 1:10:60)));
 figure
 mesh(double(h2(1:10:60, 1:10:60)));
 figure
 mesh(double(h3(1:10:60, 1:10:60)));

在这里插入图片描述

实验结论：
理想滤波器有振铃现象发生，可以看出空域滤波函数图像外围有剧烈震荡,平滑效果较为粗糙；布特沃斯滤波器几乎没有振铃现象，图像较理想滤波器而言平滑效果也比较好，高斯滤波器没有振铃现象，其平滑效果是三者最好的。

3.5高通(锐化)频域滤波器

高通滤波是低通滤波的相反过程，通过削弱傅立叶变换的低频以及保持高频相对不变来锐化图像。我们通过衰减图像的傅里叶变换的高频成分来平滑对象，因为边缘和其他灰度的急剧变化与高频分量有关，所以图像的锐化可在频域通过高通滤波来实现。一个高通滤波器是从给定的低通滤波器用下式得到：
在这里插入图片描述
其中， $H_ HP$ ( u ， v ) = 1 - $H_ LP$ ( u ， v )是低通滤波器的传递函数，也就是说之前介绍的三种低通滤波器：理想、布特沃斯、高斯低通滤波器可以通过上面的式子快速的转变为高通滤波器。

3.5.1 高通滤波函数
使用函数 lpfilter 构建另一个函数hpfilter，使用它可以产生高通滤波器:

function H = hpfilter(type,M,N,D0,n)
if nargin == 4 
    n = 1;
end
Hlp = lpfilter(type.M,N,D0,n);
H = 1 - Hlp;

我们可以观察三个滤波器的图像：

>> H = fftshift(hpfilter('ideal', 500, 500, 50)); 
>> B = fftshift(hpfilter('btw', 500, 500, 50));
>> G = fftshift(hpfilter('gaussian', 500, 500, 50));
>> 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('高斯高通滤波器图像');

在这里插入图片描述
可以从图中看出，理想高通滤波器在圆的周围明显的不连续，因此跟LPF一样会产生振铃现象。

理想高通滤波器：

 f = im2double(rgb2gray(imread('xl.jpg')));
>> PQ = paddedsize(size(f));
>> [U,V] = dftuv(PQ(1),PQ(2));
>> D0 = 0.05*PQ(2);
>> F = fft2(f, PQ(1), PQ(2));     %傅里叶变换
>> H = hpfilter('ideal',PQ(1),PQ(2),D0,2);
>> g = dftfilt(f,H);
>>subplot(2, 2, 1), imshow('xl.jpg'), title('原图像');
subplot(2, 2, 2), imshow(fftshift(H),[]), title('以图像显示的理想高通滤波器');
subplot(2, 2, 3), imshow(log(1+abs(fftshift(F))),[]), title('对数增强并居中的fft');
subplot(2, 2, 4), imshow(g,[]), title('理想高通滤波后图像的谱');

在这里插入图片描述
布特沃斯高通滤波器：

f = im2double(rgb2gray(imread('f.jpg')));
PQ = paddedsize(size(f));      %获取填充尺寸
DO = 0.05*PQ(1);               %设定滤波器的阈值
H = fftshift(hpfilter('btw', PQ(1),PQ(2), DO)); 
g = dftfilt(f,H);
figure;
subplot(1, 2, 1), imshow('f.jpg'), title('原图像');
subplot(1, 2, 2), imshow(g), title('理想高通滤波后结果');

在这里插入图片描述
高斯高通滤波器：

 f = im2double(rgb2gray(imread('xl.jpg')));
>> PQ = paddedsize(size(f));      %获取填充尺寸
>> [U,V] = dftuv(PQ(1),PQ(2));
>> D0 = 0.05*PQ(2);
>> F = fft2(f, PQ(1), PQ(2));     %傅里叶变换
>> H = fftshift(hpfilter('gaussian', PQ(1),PQ(2), D0));
>> g = dftfilt(f,H);
>> figure;
subplot(1, 2, 1), imshow('xl.jpg'), title('原图像');
subplot(1, 2, 2), imshow(g,[]), title('理想高通滤波后结果');

在这里插入图片描述
**实验总结：**以上的三个高通滤波器的共同点都是消除模糊，突出边缘，使低频分量得到了抑制，从而增强了图像的高频分量，实现图像的锐化。但是对比三个滤波器，可以明显的看出理想高通滤波器在处理过渡的频率时，可以从理想高通滤波器透视图看出它不是平滑的，因此会出现明显的振铃现象。而布特沃斯与高斯高通滤波器对图像的滤波效果都较好，都没有明显的振铃现象。其中的高斯高通滤波器滤波半径与阶数对滤波的效果也有影响，如下图在这里插入图片描述
右侧的滤波半径为5%，左侧为1%，图像的锐化程度略有提升。

3.5.2高频强调滤波
高通滤波器偏离了 dc的0项，因此减少了图像中平均值为 0的值。一种补偿方法是为高通滤波器加上偏移量。如果偏移量与将滤波器乘以某个大于 1 的常数结合起来，这种方法就叫做高频强调滤波，因为这个常量乘数突出了高频部分。这个乘数也增加了低频部分的幅度，但是只要偏移量与乘数相比较小，低频增强的影响就弱于高频增强的影响。高频强调滤波有如下传递函数：
在这里插入图片描述
a是偏移量，b 是乘数。示例：

f = im2double(rgb2gray(imread('xl.jpg')));    %把图像变为灰度图像
>> PQ = paddedsize(size(f)); 
DO = 0.05*PQ(1);    %，Do的值等于填充过 的图像垂直尺寸的5%
HBW = hpfilter('btw', PQ(1),PQ(2), DO, 2); 
H = 0.5 + 2*HBW;
gbw = dftfilt(f, HBW, 'fltpoint');       %函数 dftfilt 中用 fltpoint 选项以得到浮点结果
gbw = gscale(gbw); 
ghf = dftfilt(f, H, 'fltpoint'); 
ghf = gscale(ghf);
ghe = histeq(ghf, 256); 
subplot(2, 2, 1), imshow('xl.jpg'), title('原图像');
subplot(2, 2, 2), imshow(gbw,[]), title('布特沃斯滤波器滤波后的结果');
subplot(2, 2, 3), imshow(ghf,[]), title('高频强调的结果');
subplot(2, 2, 4), imshow(ghe,[]), title('高频强调后经直方图均衡化结果');

在这里插入图片描述
**实验结论：**清楚的骨结构和其他细节在其他任何三幅图像中是完全看不到的。最终的增强图像有些噪声，但这是在灰度级扩展时 X 光图像的典型现象。将高频强调滤波和直方图均衡结合起来得到的结果要好于单独使用任何一种方法。

3.6选择性滤波

3.6.1带阻和带通滤波器

function H = bandfilter(type,band,M,N,D0,W,n)
if nargin < 7
    n = 1;  %默认值
end
switch lower(type)
    case 'ideal'
        H = idealReject(D,D0,W);
    case 'btw'
        H = btwReject(D,D0,W,n);
    case 'gaussian'
        H = gaussianReject(D,D0,W);
    otherwise
        error('Unknown filter type.')
end

if strcmp(band,'pass')
    H = 1 - h;
end

%------------------------------------------------------------------%
function H = idealReject(D,D0,W)
RI = D <= D0 - (W/2);
RO = D >= D0 + (W/2);
H = tofloat(RO|RI);
%------------------------------------------------------------------%
function H = btwReject(D,D0,W,n)
H = 1./(1 + ((D*W)./(D.^2 - D0^2)).^2*n));
%------------------------------------------------------------------%
function H = gaussianReject(D,D0,W)
H = 1 - exp(-((D.^2 - D0^2)./(D.*W + eps)).^2);

3.6.2陷波带阻和陷波带通滤波器
陷波滤波器是更加有用的选择性滤波器。示例：

f = imread('f.jpg ');
f=rgb2gray(f);
[M,N] = size(f);
[f,revertclass] = tofloat(f);
F = fft2(f);
S = fftshift(log(1+abs(F))); %将直流分量移到频谱中心
S = gscale(S); %频谱
subplot(221),imshow(f),title(‘原图像’);
subplot(222),imshow(S),title(‘显示频谱图像’);
C1 = [99 154;128 163]; %交互式的尖峰
H1 = cnotch(‘gaussian’, ‘reject’,M,N,C1,5); %陷波滤波器
P1 = gscale(fftshift(H1).*(tofloat(S))); %计算滤波器变换的频谱
subplot(223),imshow(P1), title(‘滤波变换的频谱图像’);
g1 = dftfilt(f,H1); %滤波图像
g1=revertclass(g1);
subplot(224),imshow(g1), title(‘滤波图像’);

Daos2020

关注

0
点赞
踩
18

收藏

觉得还不错? 一键收藏
0
评论
数字图像处理第三章

二维离散傅里叶变换(一)理论推导（二）在 MATLAB 中计算及观察二维 DFT背景知识：傅里叶变换在时域和频域上都呈离散的形式，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作经过周期延拓成为周期信号再作变换。(一)理论推导令f(x,y)代表一幅大小为 MxN的数字图像，其中x=0,1,2,…,M-1, y=0,1,2,…,N-1, 由 F(u,v)表示的f(x,y)的二维离散傅立变换(DFT)可以
复制链接

扫一扫