python3-scipy-fftpack

发布时间：2021-08-19 00:13:02编辑：run阅读（3129）

傅里叶变换方法的基本思想是图像可以看作二维函数f,这个函数可以表示成在二个维度上正弦和余弦的加权和。




使用DFT可将(空间域/时域)图像中的一组灰度像素值转换为一组(频域)傅里叶系数，而且它是离散的，这是因为空间和变换变量只可以使用离散连续整数的值。频域中的傅里叶系数二维数组可以通过离散傅里叶逆变换IDFT变换至空间域，这也称为傅里叶系数重建图像。




为什么需要DFT？

  将图像变换到频域可以更好地理解图像，频域中的低频对应图像中信息的平均总体水平，而高频对应于边缘，噪声和更信息的信息。DFT对图像压缩是非常有用的，尤其是对于稀疏傅里叶图像，只有少数的傅里叶系数图像需要重建图像，因此只有那些频率可以被存储，其他的可以丢弃，这会导致高压缩。




快速傅里叶变换(FFT)是一种分而治之的算法，对于一幅N x N的图像，它递归计算DFT的速度很得多,时间复杂度O(Nlog2N)

在python中，numpy和scipy库都提供了使用FFT算法计算2D DFT/IDFT的函数。




FFT的scipy.fftpack模块

基于灰度图像的FFT算法，利用scipy.fftpack模块的fft2()/ifft2()函数来计算DFT/IDFT。

模块安装: pip install scipy

from PIL import Image
import scipy.fftpack as fp
import matplotlib.pylab as pylab
import numpy as np

def signaltonoise(a, axis=0, ddof=0):
    a = np.asanyarray(a)
    m = a.mean(axis)
    sd = a.std(axis=axis, ddof=ddof)
    return np.where(sd == 0, 0, m/sd)

# 指定默认字体
pylab.rcParams['font.sans-serif'] = ['KaiTi']
# 解决保存图像是负号'-'显示为方块的问题
pylab.rcParams['axes.unicode_minus'] = False
im = np.array(Image.open(r'D:\image_processing\image_two\hha.jpg').convert('L'))
snr = signaltonoise(im)
print('SNR for the original image = ' + str(snr))
freq = fp.fft2(im)
im1 = fp.ifft2(freq).real
snr1 = signaltonoise(im1)
print('SNR1 for the original image = ' + str(snr))
assert(np.allclose(im, im1))
pylab.figure(figsize=(20, 10))
pylab.subplot(121), pylab.imshow(im, cmap='gray'), pylab.axis('off')
pylab.title('原始图像', size=25)
pylab.subplot(122), pylab.imshow(im1, cmap='gray'), pylab.axis('off')
pylab.title('重建后的图像', size=25)
pylab.show()



从内联输出的信噪比和输入与重建图像的视觉差异可以看出，重建图像丢失了一些信息。







绘制频谱图。由于傅里叶系数是复数，因为可以直观的看出其幅度。显示傅里叶变换的幅度称变换的频谱，使用fftshift()函数使直流分量在中间.

from PIL import Image
import scipy.fftpack as fp
import matplotlib.pylab as pylab
import numpy as np

def signaltonoise(a, axis=0, ddof=0):
    a = np.asanyarray(a)
    m = a.mean(axis)
    sd = a.std(axis=axis, ddof=ddof)
    return np.where(sd == 0, 0, m/sd)

# 指定默认字体
pylab.rcParams['font.sans-serif'] = ['KaiTi']
# 解决保存图像是负号'-'显示为方块的问题
pylab.rcParams['axes.unicode_minus'] = False
im = np.array(Image.open(r'D:\image_processing\image_two\hha.jpg').convert('L'))
snr = signaltonoise(im)
print('SNR for the original image = ' + str(snr))
freq = fp.fft2(im)
freq2 = fp.fftshift(freq)
pylab.figure(figsize=(10, 10)),pylab.imshow((20*np.log10(0.1+freq2)).astype(int))
pylab.title('傅里叶频谱', size=25)
pylab.show()









FFT的numpy.fft模块

用numpy.fft模块的类似功能集计算图像的DFT,实现计算DFT的幅值和相位。利用fft2()得到傅里叶系数的实分量和虚分量，然后计算幅值，频谱和相位，最后利用ifft2()重建图像.

from skimage.io import imread
import scipy.fftpack as fp
from skimage.color import rgb2gray
import matplotlib.pylab as pylab
import numpy as np

# 指定默认字体
pylab.rcParams['font.sans-serif'] = ['KaiTi']
# 解决保存图像是负号'-'显示为方块的问题
pylab.rcParams['axes.unicode_minus'] = False
im1 = rgb2gray(imread(r'D:\image_processing\image_two\aaz.jpg'))
pylab.figure(figsize=(12,10))
freq1 = fp.fft2(im1)
im1_ = fp.ifft2(freq1).real
pylab.subplot(2,2,1)
pylab.imshow(im1, cmap='gray')
pylab.title('原始图像', size=25)
pylab.axis('off')
pylab.subplot(2,2,2)
pylab.imshow(20*np.log10(0.01+np.abs(fp.fftshift(freq1))), cmap='gray')
pylab.title('FFT频谱幅度图', size=25)
pylab.axis('off')
pylab.subplot(2,2,3)
pylab.imshow(np.angle(fp.fftshift(freq1)), cmap='gray')
pylab.title('FFT相位图', size=25)
pylab.axis('off')
pylab.subplot(2,2,4)
pylab.imshow(np.clip(im1_, 0, 255), cmap='gray')
pylab.title('重建图像', size=25)
pylab.axis('off')
pylab.show()









利用一幅图像的频谱实分量和另一幅图像的频谱虚分离来看重建的输出图像是如何变得扭曲的。

from skimage.io import imread
import scipy.fftpack as fp
from skimage.color import rgb2gray
import matplotlib.pylab as pylab
import numpy as np

# 指定默认字体
pylab.rcParams['font.sans-serif'] = ['KaiTi']
# 解决保存图像是负号'-'显示为方块的问题
pylab.rcParams['axes.unicode_minus'] = False
im1 = rgb2gray(imread(r'D:\image_processing\image_two\aaz2.jpg'))
freq1 = fp.fft2(im1)
im2 = rgb2gray(imread(r'D:\image_processing\image_two\hha.jpg'))
freq2 = fp.fft2(im2)
pylab.figure(figsize=(20,15))
im1_ = fp.ifft2(np.vectorize(complex)(freq1.real, freq2.imag)).real
im2_ = fp.ifft2(np.vectorize(complex)(freq2.real, freq1.imag)).real
pylab.subplot(221), pylab.imshow(np.clip(im1_, 0, 255), cmap='gray'), pylab.axis('off')
pylab.title('重建图像(Re(F1)+Im(F2))', size=25)
pylab.subplot(222), pylab.imshow(np.clip(im2_, 0, 255), cmap='gray'), pylab.axis('off')
pylab.title('重建图像(Re(F2)+Im(F1))', size=25)
pylab.show()

关键字：