1、基础
卷积定理
折叠误差补零
当处理DFT时，图像及其变换是周期的。在周期接近函数非零部分的持续周期时，对周期函数进行卷积会引起相邻周期的串扰。
这种称为 折叠误差 的串扰可通过补零方法来避免。
假设函数 f(x,y) 和 h(x,y) 的大小分别为 A x B 和 C x D。通过对 f 和 g 补零，构造两个大小均为 P x Q 的扩展函数。
按如下方式进行选择可以避免折叠误差：
下面 paddedsize 的函数可以满足上面等式的 P 和 Q 的最小偶数值。
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) && strcmp(PARAM, 'pwr2')
m = max([AB CD]);
P = 2^nextpow2(2*m);
PQ = [P, P];
else
error('Wrong number of inputs.')
end
示例
[f, revertcalss] = tofloat(f);
sig = 10;
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 = revertcalss(gpc);
imshow(gp)
2、DFT 滤波的基本步骤
1、使用函数 tofloat 把输入图像转换为浮点图像：
[f, revertclass] = tofloat(f);
2、使用函数 paddedsize 获得填充参数：
PQ = paddedsize(size(f));
3、得到有填充图像的傅里叶变换：
F = fft2(f, PQ(1), PQ(2));
4、生成一个大小为 PQ(1) * PQ(2) 的滤波器函数。如果它是居中的，使用前要令 H = ifftshift(H)。
5、使用滤波器乘以该变换：
G = H .* F;
6、获得G的IFFT：
g = ifft2(G);
7、将左上部的矩形修剪为原始大小：
g = g(1:size(f, 1), 1:size(f, 2));
8、需要时，将滤波后的图像转换为输入图像的类：
g = revertclass(g);
- 用于频率域滤波的M函数
function g = dftfilt(f, H, classout)
[f, revertcalss] = 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 = revertcalss(g);
elseif strcmp(classout, 'fltpoint')
return
else
error('Undefined class fot the output image.')
end
3、从空间滤波器获得频域滤波器
通常当滤波器较小时，在计算上空间滤波要比频率域滤波更有效率。
当滤波器大约有 32 个或更多元素时，使用 FFT 算法的滤波处理要比空间实现快。
函数 freqz2 可以将空间滤波器转换为频率域滤波器
H = freqz2(h, R, C)
示例
h = fspecial('sobel')'
freqz2(h)
% 频率域滤波器
PQ = paddedsize(size(f));
H = freqz2(h, PQ(1), PQ(2));
H1 = ifftshift(H)
函数 dftfilt 用于处理频率域滤波
gf = dftfilt(f, H1);
4、在频率域中直接生成滤波器
- 创建用于实现频率域滤波器的网格数组
function [U, V] = dftuv(M, N)
u = single(0:(M - 1));
v = single(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);
% ----------------------- %
[U, V] = dftuv(8, 5);
DSQ = hypot(U, V)
-
低通（平滑）频率域滤波器
理想低通滤波器（ILPF）
巴特沃斯低通滤波器（BLPF）
高斯低通滤波器
示例，对图像 f 应用一个高斯低通滤波器，其中 D0 为所填充图像宽度的 5%。
[f, revertclass] = tofloat(f);
PQ = paddedsize(size(f));
[U, V] = dftuv(PQ(1), PQ(2));
D = hypot(U, V);
D0 = 0.05 * PQ(2);
F = fft2(f, PQ(1), PQ(2));
H = exp(-(D .^ 2) / (2 * (D0^2)));
g = dftfilt(f, H);
g = revertclass(g)
- 低通滤波器整合函数
function H = lpfilter(type, M, N, D0, n)
[U, V] = dftuv(M, N);
D = hypot(U, V);
switch type
case 'ideal'
H = single(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
5、绘制线框图和表面图
绘制二维函数的最简单方法是使用函数 mesh
mesh(H(1:k:end, 1:k:end))
默认绘出彩色网格图命令为
colormap([0 0 0])
关闭网格和坐标轴命令
grid off
axis off
观察点由函数 view 控制
view(az, el) %az,el 代表方位角和仰角
view(3) %默认观察点
示例，绘制一个高斯低通滤波器
H = fftshift(test2('gaussian', 500, 500, 50));
mesh(double(H(1:10:500, 1:10:500)))
axis tight
colormap([0 0 0])
- 表面图
surf(H)
6、高通（锐化）频率域滤波器
理想高通滤波器（IHPF）
巴特沃斯高通滤波器（BHPF）
高斯高通滤波器
基于 lpfilter 构建一个高通滤波器函数
function H = hpfilter(type, M, N, D0, n)
if nargin == 4
n = 1;
end
H1p = lpfilter(type, M, N, D0, n);
H = 1 - H1p;
% ---------------------------%
H = fftshift(hpfilter('ideal', 500, 500, 50));
高频强调滤波
高通滤波器偏离了直流项，将图像的平均值降低为0.补偿方法之一是给高通滤波器加上一个偏移量。把偏移量与将滤波器乘以一个大于1的常数结合起来的方法就称为 高频强调滤波。
示例，联合使用高频强调滤波和直方图均衡
PQ = paddedsize(size(f));
D0 = 0.05 * PQ(1);
HBW = hpfilter('btw', PQ(1), PQ(2), D0, 2);
H = 0.5 + 2*HBW;
gbw = dftfilt(f, HBW, 'fltpoint');
gbw = gscale(gbw);
ghf = dftfilt(f, H, 'fltpoint');
ghf = gscale(ghf);
ghe = histeq(ghf, 256);
