Python OpenCV 图像处理之傅里叶变换,取经之旅第 52 篇

举报
梦想橡皮擦 发表于 2021/09/24 14:10:42 2021/09/24
【摘要】 Python OpenCV 365 天学习计划,与橡皮擦一起进入图像领域吧。本篇博客是这个系列的第 52 篇。 学在前面傅里叶变换(Fourier Transform,FT)今天第一次接触,按照以往的套路,我们依旧是先应用起来,然后逐步的迭代学习。傅里叶变换是对同一个事物观看角度的变化,不在从时域进行观看,从频域进行观看。这里提及了两个新词,时域和频域,先简单理解一下,时域,在时间范围内事物...

Python OpenCV 365 天学习计划,与橡皮擦一起进入图像领域吧。本篇博客是这个系列的第 52 篇。

学在前面

傅里叶变换(Fourier Transform,FT)今天第一次接触,按照以往的套路,我们依旧是先应用起来,然后逐步的迭代学习。

傅里叶变换是对同一个事物观看角度的变化,不在从时域进行观看,从频域进行观看。这里提及了两个新词,时域和频域,先简单理解一下,时域,在时间范围内事物发生的变化,频域是指的事物的变化频率,是不是很绕,没啥问题,先用,先会调 API。

傅里叶原理略过,先说一下应用它之后对图像产生的影响。

  • 使用高通滤波器之后,会保留高频信息,增强图像细节,例如边界增强;
  • 使用低通滤波器之后,会保留低频信息,边界模糊。

傅里叶变换应用

原理既然先放下了,那先应用起来吧,我们将分别学习 numpy 和 OpenCV 两种方式实现傅里叶变换。

Numpy 实现傅里叶变换

通过 numpy 实现傅里叶变换用到的函数是 np.fft.fft2 ,函数原型如下:

fft2(a, s=None, axes=(-2, -1), norm=None)

参数说明如下:

  • a:输入图像;
  • s:整数序列,输出数组的大小;
  • axex:整数序列,用于计算 FFT 的可选轴;
  • norm:规范化模式。

有些参数如果不实际使用,比较难看出结果,当然函数的说明,还是 官方手册 最靠谱。应用层面的代码编写流程如下:
通过 np.fft.fft2 函数进行傅里叶变换得到频率分布,再调用 np.fft.fftshift 函数将中心位置转移至中间。

测试代码如下:

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

img = cv.imread('test.jpg', 0)

# 快速傅里叶变换算法得到频率分布,将空间域转化为频率域
f = np.fft.fft2(img)

# 默认结果中心点位置是在左上角,通过下述代码将中心点转移到中间位置
# 将低频部分移动到图像中心
fshift = np.fft.fftshift(f)

# fft 结果是复数, 其绝对值结果是振幅
result = 20*np.log(np.abs(fshift))

plt.subplot(121)
plt.imshow(img, cmap='gray')
plt.title('original')
plt.axis('off')

plt.subplot(122)
plt.imshow(result, cmap='gray')
plt.title('result')
plt.axis('off')

plt.show()

这个是很多地方反复提及的案例,在补充一下相关内容。
np.fft.fft2 是一个频率转换函数,它的第一个参数是输入图像,即灰度图像。第二个参数是可选的,它决定输出数组的大小。
第二个参数的大小如果大于输入图像的大小,则在计算 FFT 之前用零填充输入图像;如果小于输入图像,将裁切输入图像。如果未传递任何参数,则输出数组的大小将与输入的大小相同,测试如下:

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

# 生成一个图片
src = np.ones((5, 5), dtype=np.uint8)*100

print(src)
print(src.shape)

f = np.fft.fft2(src,(7,7))
print(f.shape)
print(f)

中心点进行移动之后,可以参考下图运行结果,而且可以知道 np.fft.fft2 函数应用之后得到的是复数矩阵。

后面要进行的操作,有部分数学知识,暂未搞懂,摘录如下:
进行完频率变换之后,就可以构建振幅谱了,最后在通过对数变换来压缩范围。
或者可以理解为将复数转为浮点数进行傅里叶频谱图显示,补充代码如下:

fimg = np.log(np.abs(fshift))

最终得到的结果如下,当然这个是采用的随机图,如果换成一张灰度图,可以验证如下。

修改之后的代码如下:

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

# 生成一个图片
# src = np.ones((5, 5), dtype=np.uint8)*100
src = cv.imread("./test.jpg", 0)
print("*"*100)
print(src)
print(src.shape)

f = np.fft.fft2(src)
print("*"*100)
print(f)

fshift = np.fft.fftshift(f)
print("*"*100)
print(fshift)
# 将复数转为浮点数进行傅里叶频谱图显示
fimg = 20*np.log(np.abs(fshift))
print(fimg)

# 图像显示
plt.subplot(121), plt.imshow(src, "gray"), plt.title('origin')
plt.axis('off')
plt.subplot(122), plt.imshow(fimg, "gray"), plt.title('Fourier')
plt.axis('off')
plt.show()

基本结论:左边为原始灰度图像,右边为频率分布图谱,其越靠近中心位置频率越低,灰度值越高亮度越亮的中心位置代表该频率的信号振幅越大。
接下来将其进行逆变换,也就是反向操作回去,通过 numpy 实现傅里叶逆变换,它是傅里叶变换的逆操作,将频谱图像转换为原始图像。用到核心函数与原型分别如下。
np.fft.ifft2

# 实现图像逆傅里叶变换,返回一个复数数组
np.fft.ifft2(a, s=None, axes=(-2, -1), norm=None)

np.fft.ifftshift

# fftshit()函数的逆函数,它将频谱图像的中心低频部分移动至左上角
np.fft.ifftshift(x, axes=None)

依据上述内容,对傅里叶变换进行逆向的操作,测试代码如下:

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

# 生成一个图片
# src = np.ones((5, 5), dtype=np.uint8)*100
src = cv.imread("./test.jpg", 0)
print("*"*100)
print(src)
print(src.shape)

f = np.fft.fft2(src)
print("*"*100)
print(f)

fshift = np.fft.fftshift(f)
print("*"*100)
print(fshift)

# 将复数转为浮点数进行傅里叶频谱图显示
fimg = np.log(np.abs(fshift))
print(fimg)

# 逆傅里叶变换
ifshift = np.fft.ifftshift(fshift)
# 将复数转为浮点数进行傅里叶频谱图显示
ifimg = np.log(np.abs(ifshift))
if_img = np.fft.ifft2(ifshift)

origin_img = np.abs(if_img)

# 图像显示
plt.subplot(221), plt.imshow(src, "gray"), plt.title('origin')
plt.axis('off')
plt.subplot(222), plt.imshow(fimg, "gray"), plt.title('fourier_img')
plt.axis('off')
plt.subplot(223), plt.imshow(origin_img, "gray"), plt.title('origin_img')
plt.axis('off')
plt.subplot(224), plt.imshow(ifimg, "gray"), plt.title('ifimg')
plt.axis('off')
plt.show()

最终的结果如下:

如果在上述傅里叶变换之后的频谱图像中,增加一个低通滤波,将会得到如下结果。

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

# 实现傅里叶变换的低通滤波
src = cv.imread("./test.jpg", 0)
f = np.fft.fft2(src)
fshift = np.fft.fftshift(f)

rows, cols = src.shape
crow, ccol = rows//2, cols//2

# 制定掩模,大小和图像一样,np.zeros初始化
mask = np.zeros((rows, cols), np.uint8)
# 使中心位置,上下左右距离30,置为1
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# 掩模与 DFT 后结果相乘只保留出中间区域
fshift = fshift*mask

# 逆傅里叶变换
ifshift = np.fft.ifftshift(fshift)
# 将复数转为浮点数进行傅里叶频谱图显示
ifimg = np.fft.ifft2(ifshift)
dft_img = np.abs(ifimg)


# 图像显示
plt.subplot(121), plt.imshow(src, "gray"), plt.title('origin')
plt.axis('off')
plt.subplot(122), plt.imshow(dft_img, "gray"), plt.title('dft_img')
plt.axis('off')
plt.show()

如果希望实现高通滤波,只需要修改掩模数据即可。

mask = np.ones((rows, cols), np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 0

参考了很多文章,但是一篇文章被反复提及,这里也备注一下:https://zhuanlan.zhihu.com/p/19763358
当然傅里叶变换在官方手册的内容也需要时长进行阅读一下下,官网地址

OpenCV 实现傅里叶变换

OpenCV 实现傅里叶变换与 numpy 基本一致,核心差异在函数的使用上,具体区别如下:
Opencv 中主要通过 cv2.difcv2.idif(傅里叶逆变换),在输入图像之前需要先转换成把图像从 np.uint8 转换为 np.float32 格式,其得到的结果中频率为 0 的部分会在左上角,要转换到中心位置,通过 shift 变换来实现,与 numpy 一致,cv2.dif 返回的结果是双通道的(实部,虚部),还需要转换成图像格式才能展示。

函数 cv2.dft() 原型如下

dst = cv.dft(src[, dst[, flags[, nonzeroRows]]])

参数说明:

  • src:输入图像,要求 np.float32 格式;
  • dst:输出图像,双通道(实部+虚部),大小和类型取决于第三个参数 flags
  • flags:表示转换标记,默认为 0,存在多种取值,参见后文;
  • nonzeroRows:默认为 0,暂时不考虑。

flags 取值如下:

  • DFT_INVERSE:用一维或二维逆变换取代默认的正向变换;
  • DFT_SCALE: 缩放比例标识符,根据数据元素个数平均求出其缩放结果,如有 N 个元素,则输出结果以 1/N 缩放输出,常与 DFT_INVERSE 搭配使用;
  • DFT_ROWS:对输入矩阵的每行进行正向或反向的傅里叶变换;此标识符可在处理多种适量的的时候用于减小资源的开销,这些处理常常是三维或高维变换等复杂操作;
  • DFT_COMPLEX_OUTPUT:对一维或二维的实数数组进行正向变换,这样的结果虽然是复数阵列,但拥有复数的共轭对称性(CCS),可以以一个和原数组尺寸大小相同的实数数组进行填充,这是最快的选择也是函数默认的方法。你可能想要得到一个全尺寸的复数数组(像简单光谱分析等等),通过设置标志位可以使函数生成一个全尺寸的复数输出数组;
  • DFT_REAL_OUTPUT:对一维二维复数数组进行逆向变换,这样的结果通常是一个尺寸相同的复数矩阵,但是如果输入矩阵有复数的共轭对称性(比如是一个带有 DFT_COMPLEX_OUTPUT标识符的正变换结果),便会输出实数矩阵。

以上内容摘抄网络,原创人已经无法找到,尴尬。

总结下来就是:

  • DFT_COMPLEX_OUTPUT:得到一个复数形式的矩阵;
  • DFT_REAL_OUTPUT:只输出复数的实部;
  • DFT_INVERSE:进行傅里叶逆变换;
  • DFT_SCALE:是否除以 MxN (M 行 N 列的图片,共有有 MxN 个像素点);
  • DFT_ROWS:输入矩阵的每一行进行傅里叶变换或者逆变换。

最后需要注意的是,输出的频谱结果是一个复数,需要调用 cv2.magnitude() 函数将傅里叶变换的双通道结果转换为 0 到 255 的范围。

这个函数比较简单,原型和说明如下。

cv2.magnitude 函数原型:cv2.magnitude(x, y)

  • x 表示浮点型 X 坐标值,即实部
  • y 表示浮点型 Y 坐标值,即虚部

基础知识简单过一遍,直接进行案例的说明。

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']

src = cv.imread("test.jpg", 0)

# OpneCV傅里叶变换函数
# 需要将图像进行一次float转换
result = cv.dft(np.float32(src), flags=cv.DFT_COMPLEX_OUTPUT)
# 将频谱低频从左上角移动至中心位置
dft_shift = np.fft.fftshift(result)
# 频谱图像双通道复数转换为 0-255 区间
result1 = 20 * np.log(cv.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
# 图像显示
plt.subplot(121), plt.imshow(src, 'gray'), plt.title('原图像')
plt.axis('off')
plt.subplot(122), plt.imshow(result1, 'gray'), plt.title('傅里叶变换')
plt.axis('off')
plt.show()

运行结果与 numpy 一致。

套用上面的知识,使用 OpenCV 里面的函数对图片使用低通滤波。

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']

src = cv.imread("test.jpg", 0)

# OpneCV傅里叶变换函数
# 需要将图像进行一次float转换
result = cv.dft(np.float32(src), flags=cv.DFT_COMPLEX_OUTPUT)
# 将频谱低频从左上角移动至中心位置
dft_shift = np.fft.fftshift(result)
# 频谱图像双通道复数转换为 0-255 区间
result = 20 * np.log(cv.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))

rows, cols = src.shape
crow, ccol = rows//2, cols//2

mask = np.zeros((rows, cols, 2), np.uint8)
mask[int(crow-30):int(crow+30), int(ccol-30):int(ccol+30)] = 1
#  LPF(低通滤波)
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv.idft(f_ishift)
img_back = cv.magnitude(img_back[:, :, 0], img_back[:, :, 1])

# 图像显示
plt.subplot(131), plt.imshow(src, 'gray'), plt.title('原图像')
plt.axis('off')
plt.subplot(132), plt.imshow(result, 'gray'), plt.title('傅里叶变换')
plt.axis('off')
plt.subplot(133), plt.imshow(img_back, 'gray'), plt.title('低通滤波之后的图像')
plt.axis('off')
plt.show()

高通滤波的代码和运行效果,就交给你自己来实现吧。

橡皮擦的小节

HPF(高通滤波),LPF(低通滤波)
希望今天的 1 个小时(貌似不太够)你有所收获,我们下篇博客见~

相关阅读


技术专栏

  1. Python 爬虫 100 例教程,超棒的爬虫教程,立即订阅吧
  2. Python 爬虫小课,精彩 9 讲

今天是持续写作的第 97 / 100 天。
如果你想跟博主建立亲密关系,可以关注同名公众号 梦想橡皮擦,近距离接触一个逗趣的互联网高级网虫。
博主 ID:梦想橡皮擦,希望大家点赞评论收藏

【版权声明】本文为华为云社区用户原创内容,未经允许不得转载,如需转载请自行联系原作者进行授权。如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱: cloudbbs@huaweicloud.com
  • 点赞
  • 收藏
  • 关注作者

评论(0

0/1000
抱歉,系统识别当前为高风险访问,暂不支持该操作

全部回复

上滑加载中

设置昵称

在此一键设置昵称,即可参与社区互动!

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。