注意:这一系列实验的图像处理程序,使用Matlab实现最重要的图像处理算法
1.Fourier兑换
(1)频域增强
除了在空间域内能够加工处理图像以外,我们还能够将图像变换到其它空间后进行处理。这些方法称为变换域方法,最常见的变换域是频域。
使用Fourier变换把图像从空间域变换到频域。在频域内做对应增强处理,再从频域变换到空间域得到处理后的图像。
我们这里主要学习Fourier变换和FFT变换的算法,没有学过通信原理,我对信号、时域分析也不是非常清楚。
2.FFT算法
(1)离散Fourier变换,DFT
函数
非常显然求出全部的长度为N的信号的DFT须要
(2)高速Fourier变换,FFT
高速傅立叶变换FFT是利用单位复数根的特殊性质(消去引理和折半引理,见《算法导论》(第3版中文版)P532具体论述)。在
FFT有两种基本实现方式:
- 递归FFT
- 迭代FFT。也叫FFT蝶式运算
递归FFT因为递归栈开销大且容量有限等缺陷(但理解easy),一般计算採取迭代FFT实现。
(3)迭代FFT
直接给出算法导论版本号的迭代FFT算法:
当中求逆序数拷贝的函数为:
FFT採用折半迭代的思想,因此速度能减少到
(4)迭代FFTMatlab实现
Matlab有fft函数,我们也能够自己实现:
function [ fft_m ] = IterativeFFT( vec )
clear i;
n = length(vec);
fft_m = BitReverseCopy(vec);
for s = 1 : log2(n)
m = power(2, s);
wm = exp(- 2 * pi * i / m);
for k = 0 : m : n - 1
w = 1;
for j = 0 : m / 2 - 1
t = w * fft_m(k + j + m / 2 + 1);
u = fft_m(k + j + 1);
fft_m(k + j + 1) = u + t;
fft_m(k + j + m / 2 + 1) = u - t;
w = w * wm;
end
end
end
end
BitReverseCopy函数例如以下:
function [ copy ] = BitReverseCopy( vec )
n = length(vec);
copy = zeros(1, n);
bitn = log2(n);
for i = 0 : n - 1
revi = bin2dec(fliplr(dec2bin(i, bitn)));
copy(revi + 1) = vec(i + 1);
end
end
须要特别注意的是:
- 一般给出的FFT算法伪代码都是基于下标从零開始的数组。写在Matlab须要考虑映射关系
clear i
是为了怕之前有变量i和复数记号i混淆,清楚Matlab workspace中的缓存- 默认vec是double类型的。因为中间採用很多double类型运算
3.图像的二维Fourier变换
(1)二维DFT
二维DFT定义公式为:
计算一个频域点须要
(2)将二维DFT分解为两个一维DFT
Fourier变换的变换核(公式中和
先对二维矩阵的每一列做一维DFT:
再对变换后的矩阵的每一行做一维DFT:
最后以
(3)用一维FFT实现二维FFT
相同的我们能够用两个一维FFT实现二维FFT,时间复杂度
function [ mfft2 ] = JCGuoFFT2( data )
h = size(data, 1);
w = size(data, 2);
mfft2 = data;
if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w
disp(\'JCGuoFFT2 exit: h and w must be the power of 2!\')
else
for i = 1 : h
mfft2(i, :) = IterativeFFT(mfft2(i, :));
end
for j = 1 : w
mfft2(:, j) = IterativeFFT(mfft2(:, j));
end
end
end
代码非常简单,先做FFT行变换再做FFT列变换。之前忘记提到。我这里实现的都是长度必须是2的次方的Fourier变换。因此有时候会做一些长度规范检查。
(4)变换结果
经过JCGuoFFT2的二维傅立叶变换,我们能够得到复平面内各个点的复数值,那么显示在图像中须要先求出幅值:
pic1_fft = JCGuoFFT2(double(pic1));
pic1_fft_amp = abs(pic1_fft);
在对幅值做一次log变换得到较好的频域图像:
pic1_fft_amp_log = log(1 + pic1_fft_amp);
绘制结果例如以下图:
(5)低频信号移到图像中心点
因为复数运算的周期特性,图像的Fourier变换在复平面上是全然对称的。能够想象平面是由无限多个上图(右)频域图像拼接而成的二维阵列。一般研究频域图像是把低频部分,也就是变换后的边缘部分移到图像中心点。Matlab提供fftshift
函数完毕平移。
平移的思路有两个
- 通过Fourier变换平移定理先把原始图像做变换再做FFT
- 先做FFT后再依据频域图像的对称性做对称变换
查阅Matlab文档发现它是採用另外一种方法,对图像做下面子矩阵交换:
那么我们能够非常easy的写出自己的fftshift
:
function [ after ] = FFTShift( before )
h = size(before, 1);
w = size(before, 2);
after = before;
if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w
disp(\'FFTShift exit: h and w must be the power of 2!\')
else
hh = h / 2;
hw = w / 2;
after(1 : hh, 1 : hw) = before(hh + 1 : h, hw + 1 : w);
after(1 : hh, hw + 1 : w) = before(hh + 1 : h, 1 : hw);
after(hh + 1 : h,
全部评论
请发表评论