天津中医药大学-20级医学信息工程 教师：王翌 学生：邓集亲
声明：本文章中所涉及的代码并非我个人独立完成的成果。
在撰写的过程中，我广泛地吸取了前一辈人，尤其是学长学姐们的宝贵学习经验。通过深入研究他们的学习轨迹，以及查阅和分析了众多相关的理论文献与资料，并在王老师的悉心指导下，经过反复的实验、调试与优化，最终得以总结完成本文所展现的代码。
建议各位学弟学妹们，先根据王老师的授课内容，独立思考一下应该如何完成。如果实在是有理解上的困难，不知道从何下手，再来学习和参考本文。
学长我是用Python写的，如果你使用MATLAB，同样可以参考此代码进行翻译。我在代码中加入了许多注释和测试环节，以便于理解。由于学长能力有限，代码中或许存在一些疏漏或错误，请批判性地参考。
实验四：FFT变换
作业要求：
- 参考“傅里叶变换”课的内容，对输入图像进行快速傅里叶变换。
- 通过傅里叶逆变换，还原出原图片。
实验图片：
fft.tif，再自选其它至少20张图片（包括灰度图片和彩色图片）。
傅里叶变换用于分析图像的频率特性，对采样的离散信号进行傅里叶变换，得到周期性的N点DFT。可以将图像视为在两个方向上采样的信号。因此，在X和Y方向都进行傅立叶变换，可以得到图像的频率表示。
更直观地说，对于正弦信号，如果幅度在短时间内变化很快，则可以说它是高频信号。如果变化缓慢，则为低频信号。
将相同的想法扩展到图像。图像中的振幅在哪里急剧变化？在边缘点或噪声。因此，可以说边缘和噪声是图像中的高频内容。如果幅度没有太大变化，则它是低频分量。
在图像处理中，频域反应了图像在空域灰度变化剧烈程度，也就是图像灰度的变化速度，也就是图像的梯度大小。也就是说，傅立叶变换提供另外一个角度来观察图像， 可以将图像从灰度分布转化到频率分布上来观察图像的特征。书面一点说就是，傅里叶变换提供了一条从空域到频率自由转换的途径。
经过傅里叶变换后的图像，我们可以得到一张频幅图，图像信息只有频率、幅度、方向。
同理，逆变换则可以根据频幅图的信息还原原图像。
Windows图片文件目录：
ImageSet文件目录下储存其他待处理图片
Python代码：
#导入库
import cv2 as cv
import numpy as np
from matplotlib import pyplot
from numpy.core import integer, asarray, roll,zeros, swapaxes, take
from numpy.core.multiarray import normalize_axis_index
from numpy.fft import _pocketfft_internal as pfi
#定义函数进行fft变换计算
def raw_fft(a, n, axis, is_real, is_forward, inv_norm): #为了变换结果，这里需要定义浮点数inv_norm
axis = normalize_axis_index(axis, a.ndim)
if n is None:
n = a.shape[axis]
fct = 1/inv_norm
if a.shape[axis] != n:
s = list(a.shape)
index = [slice(None)]*len(s)
if s[axis] > n:
index[axis] = slice(0, n)
a = a[tuple(index)]
else:
index[axis] = slice(0, s[axis])
s[axis] = n
z = zeros(s, a.dtype.char)
z[tuple(index)] = a
a = z
if axis == a.ndim-1:
r = pfi.execute(a, is_real, is_forward, fct)
else:
a = swapaxes(a, axis, -1)
r = pfi.execute(a, is_real, is_forward, fct)
r = swapaxes(r, axis, -1)
return r
def fft_shift(x, axes=None): #将图像低频部分转移到图像的中心
x = asarray(x)
if axes is None:
axes = tuple(range(x.ndim))
shift = [dim // 2 for dim in x.shape]
elif isinstance(axes, integer):
shift = x.shape[axes] // 2
else:
shift = [x.shape[ax] // 2 for ax in axes]
return roll(x, shift, axes)
def ifft_shift(x, axes=None): #将图像低频部分移到原来位置
x = asarray(x)
if axes is None:
axes = tuple(range(x.ndim))
shift = [-(dim // 2) for dim in x.shape]
elif isinstance(axes, integer):
shift = -(x.shape[axes] // 2)
else:
shift = [-(x.shape[ax] // 2) for ax in axes]
return roll(x, shift, axes)
def fft(a, n=None, axis=-1): #定义一维傅里叶变换函数
a = asarray(a)
if n is None:
n = a.shape[axis]
inv_norm = 1
output = raw_fft(a, n, axis, False, True, inv_norm)
return output
def ifft(a, n=None, axis=-1):
a = asarray(a)
if n is None:
n = a.shape[axis]
inv_norm = n
output = raw_fft(a, n, axis, False, False, inv_norm)
return output
def cook_nd_args(a, shape=None, axes=None, invreal=0):
if shape is None:
shapeless = 1
if axes is None:
shape = list(a.shape)
else:
shape = take(a.shape, axes)
else:
shapeless = 0
shape = list(shape)
if axes is None:
axes = list(range(-len(shape), 0))
if invreal and shapeless:
shape[-1] = (a.shape[axes[-1]] - 1) * 2
return shape, axes
def raw_fftnd(a, shape=None, axes=None, function=fft):
a = asarray(a)
shape, axes = cook_nd_args(a, shape, axes)
itl = list(range(len(axes)))
itl.reverse()
for ii in itl:
a = function(a, n=shape[ii], axis=axes[ii])
return a
def imagefft2(a, shape=None, axes=(-2, -1)): #图片二维傅里叶变换
return raw_fftnd(a, shape, axes, fft)
def imageifft2(a, shape=None, axes=(-2, -1)):
return raw_fftnd(a, shape, axes, ifft)
#读取图像
original = cv.imread(r'D:\deng\smalljob\ImageSet\fft.tif')
original_gray = cv.cvtColor(original, cv.COLOR_BGR2GRAY) #将图像从RGB颜色空间转换到灰度颜色空间
#傅里叶变换
f = imagefft2(original_gray)
fshift = fft_shift(f) #将变换后图像的低频部分转移到图像的中心
result = np.log(np.abs(fshift)+1e-5) #改变浮点数的精度为1e-5
#傅里叶逆变换
ishift = ifft_shift(fshift)
iImg = imageifft2(ishift)
iImg = np.abs(iImg)
#展示结果
pyplot.subplot(131),pyplot.imshow(original_gray,cmap='gray')
pyplot.title('original')
pyplot.axis('off')
pyplot.subplot(132),pyplot.imshow(result,cmap='gray')
pyplot.title('fft')
pyplot.axis('off')
pyplot.subplot(133),pyplot.imshow(iImg,cmap='gray')
pyplot.title('ifft')
pyplot.axis('off')
pyplot.show()
运行结果：
以下是对彩色图像进行傅里叶变换，即拆分图像颜色通道，分别进行傅里叶变换
#导入库
import cv2 as cv
import numpy as np
from matplotlib import pyplot
from numpy.core import integer, asarray, roll,zeros, swapaxes, take
from numpy.core.multiarray import normalize_axis_index
from numpy.fft import _pocketfft_internal as pfi
#定义函数进行fft变换计算
def raw_fft(a, n, axis, is_real, is_forward, inv_norm): #为了变换结果，这里需要定义浮点数inv_norm
axis = normalize_axis_index(axis, a.ndim)
if n is None:
n = a.shape[axis]
fct = 1/inv_norm
if a.shape[axis] != n:
s = list(a.shape)
index = [slice(None)]*len(s)
if s[axis] > n:
index[axis] = slice(0, n)
a = a[tuple(index)]
else:
index[axis] = slice(0, s[axis])
s[axis] = n
z = zeros(s, a.dtype.char)
z[tuple(index)] = a
a = z
if axis == a.ndim-1:
r = pfi.execute(a, is_real, is_forward, fct)
else:
a = swapaxes(a, axis, -1)
r = pfi.execute(a, is_real, is_forward, fct)
r = swapaxes(r, axis, -1)
return r
def fft_shift(x, axes=None): #将图像低频部分转移到图像的中心
x = asarray(x)
if axes is None:
axes = tuple(range(x.ndim))
shift = [dim // 2 for dim in x.shape]
elif isinstance(axes, integer):
shift = x.shape[axes] // 2
else:
shift = [x.shape[ax] // 2 for ax in axes]
return roll(x, shift, axes)
def ifft_shift(x, axes=None): #将图像低频部分移到原来位置
x = asarray(x)
if axes is None:
axes = tuple(range(x.ndim))
shift = [-(dim // 2) for dim in x.shape]
elif isinstance(axes, integer):
shift = -(x.shape[axes] // 2)
else:
shift = [-(x.shape[ax] // 2) for ax in axes]
return roll(x, shift, axes)
def fft(a, n=None, axis=-1): #定义一维傅里叶变换函数
a = asarray(a)
if n is None:
n = a.shape[axis]
inv_norm = 1
output = raw_fft(a, n, axis, False, True, inv_norm)
return output
def ifft(a, n=None, axis=-1):
a = asarray(a)
if n is None:
n = a.shape[axis]
inv_norm = n
output = raw_fft(a, n, axis, False, False, inv_norm)
return output
def cook_nd_args(a, shape=None, axes=None, invreal=0):
if shape is None:
shapeless = 1
if axes is None:
shape = list(a.shape)
else:
shape = take(a.shape, axes)
else:
shapeless = 0
shape = list(shape)
if axes is None:
axes = list(range(-len(shape), 0))
if invreal and shapeless:
shape[-1] = (a.shape[axes[-1]] - 1) * 2
return shape, axes
def raw_fftnd(a, shape=None, axes=None, function=fft):
a = asarray(a)
shape, axes = cook_nd_args(a, shape, axes)
itl = list(range(len(axes)))
itl.reverse()
for ii in itl:
a = function(a, n=shape[ii], axis=axes[ii])
return a
def imagefft2(a, shape=None, axes=(-2, -1)): #图片二维傅里叶变换
return raw_fftnd(a, shape, axes, fft)
def imageifft2(a, shape=None, axes=(-2, -1)):
return raw_fftnd(a, shape, axes, ifft)
#读取图像
original = cv.imread(r'D:\deng\smalljob\ImageSet\yz.jpg')
original = cv.cvtColor(original, cv.COLOR_BGR2RGB)
b,g,r=cv.split(original)
b_1=imagefft2(b)
g_1=imagefft2(g)
r_1=imagefft2(r)
b_fshift=fft_shift(b_1)
g_fshift=fft_shift(g_1)
r_fshift=fft_shift(r_1)
result_1 = np.log(np.abs(b_fshift)+1e-5)
result_2 = np.log(np.abs(g_fshift)+1e-5)
result_3 = np.log(np.abs(r_fshift)+1e-5)
b_ishift = ifft_shift(b_fshift)
g_ishift = ifft_shift(g_fshift)
r_ishift = ifft_shift(r_fshift)
b_2 = imageifft2(b_ishift)
g_2 = imageifft2(g_ishift)
r_2 = imageifft2(r_ishift)
img=np.dstack((b_2,g_2,r_2))
img=np.abs(img).astype(np.float32)/255.0
#展示结果
pyplot.subplot(231),pyplot.imshow(original,cmap='gray')
pyplot.title('original')
pyplot.axis('off')
pyplot.subplot(232),pyplot.imshow(result_1,cmap='gray')
pyplot.title('b_fft')
pyplot.axis('off')
pyplot.subplot(233),pyplot.imshow(result_2,cmap='gray')
pyplot.title('g_fft')
pyplot.axis('off')
pyplot.subplot(234),pyplot.imshow(result_3,cmap='gray')
pyplot.title('r_fft')
pyplot.axis('off')
pyplot.subplot(235),pyplot.imshow(img,cmap='gray')
pyplot.title('ifft')
pyplot.axis('off')
pyplot.show()
运行结果：
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
