python3-scipy低通滤波器

发布时间:2021-09-01 23:30:01编辑:run阅读(5002)

    python3-scipy低通滤波器只允许图像(通过DFT获得)的频域表示的低频分量通过,并阻止超过截止值的全部高频。利用离散傅里叶逆变换重建图像,由于高频分量对应于边缘,细节,噪声等,低通滤波器往往会滤除它们。利用numpy,scipy和scikit-image库的不同函数来实现低通滤波器以及观察低通滤波器对图像的影响。


    使用Scipy ndimage和numpy fft实现LPF.

    numpy fft模块的fft2()函数也可用于图像的FFT。Scipy的ndimage模块提供了一组用于在频域中对图像应用低通滤波器的函数。fourier_gaussian()函数可实现多维高斯函数傅里叶滤波,将频率数组与已知大小的高斯核的傅里叶变换相乘。


    低通滤波器(加权均值滤波器)来模糊灰度图像.

    from skimage.io import imread
    import scipy.fftpack as fp
    import matplotlib.pylab as pylab
    import numpy as np
    from scipy import signal, misc, ndimage
    
    pylab.rcParams['font.sans-serif'] = ['KaiTi']
    # 解决保存图像是负号'-'显示为方块的问题
    pylab.rcParams['axes.unicode_minus'] = False
    fig, (axes1, axes2) = pylab.subplots(1, 2, figsize=(20, 10))
    pylab.gray()
    im = np.mean(imread(r'D:\image_processing\image3\w.jpg'), axis=2)
    freq = fp.fft2(im)
    freq_gaussian = ndimage.fourier_gaussian(freq, sigma=4)
    im1 = fp.ifft2(freq_gaussian)
    axes1.imshow(im)
    axes1.set_title('原图灰度图像', size=20)
    axes1.axis('off')
    axes2.imshow(im1.real)
    axes2.set_title('低通滤波器模糊后图像', size=20)
    axes2.axis('off')
    pylab.show()

    image.png

    显示高斯核图像频谱

    from skimage.io import imread
    import scipy.fftpack as fp
    import matplotlib.pylab as pylab
    import numpy as np
    import numpy.fft
    from scipy import signal, misc, ndimage
    
    pylab.rcParams['font.sans-serif'] = ['KaiTi']
    # 解决保存图像是负号'-'显示为方块的问题
    pylab.rcParams['axes.unicode_minus'] = False
    fig, (axes1, axes2) = pylab.subplots(1, 2, figsize=(20, 10))
    pylab.gray()
    im = np.mean(imread(r'D:\image_processing\image3\w.jpg'), axis=2)
    freq = fp.fft2(im)
    freq_gaussian = ndimage.fourier_gaussian(freq, sigma=4)
    im1 = fp.ifft2(freq_gaussian)
    
    pylab.figure(figsize=(10, 10))
    pylab.imshow((20*np.log10(0.1+numpy.fft.fftshift(freq_gaussian))).astype(int))
    pylab.title('高斯核频谱', size=25)
    pylab.show()

    image.png



    使用scipy fftpack实现低通滤波器.

    1 利用scipy.fftpack fft2实现二维快速傅里叶变换,并获得图像的频域表示。

    2 只保留低频分量(去除高频分量)

    3 执行快速傅里叶变换,以重建图像

    高频分量更多的是对应于图像的平均信息,而随着我们去除越来越多的高频分量,图像的细节信息(例如边缘)就会丢失。

    from PIL import Image
    import scipy.fftpack as fp
    from scipy import stats, fftpack
    import matplotlib.pylab as pylab
    import numpy as np
    
    # 指定默认字体
    pylab.rcParams['font.sans-serif'] = ['KaiTi']
    # 解决保存图像是负号'-'显示为方块的问题
    pylab.rcParams['axes.unicode_minus'] = False
    im = np.array(Image.open(r'D:\image_processing\image3\q1.jpg').convert('L'))
    freq = fp.fft2(im)
    (w, h) = freq.shape
    half_w, half_h = int(w/2), int(h/2)
    freq1 = np.copy(freq)
    freq2 = fftpack.fftshift(freq1)
    freq2_low = np.copy(freq2)
    freq2_low[half_w-10:half_w+11, half_h-10:half_h+11] = 0
    freq2 -= freq2_low
    im1 = fp.ifft2(fftpack.ifftshift(freq2)).real
    pylab.figure(figsize=(20,15))
    pylab.subplot(131)
    pylab.imshow(im)
    pylab.title('原始灰度图像', size=25)
    pylab.axis('off')
    pylab.subplot(132)
    pylab.imshow(im1, cmap='gray')
    pylab.title('LPF后的图像', size=25)
    pylab.axis('off')
    pylab.show()

    image.png


    在阻断高频后如何在对数域中绘制图像的频谱,换言之,就是只允许低频通过。

    from PIL import Image
    import scipy.fftpack as fp
    from scipy import stats, fftpack
    import matplotlib.pylab as pylab
    import numpy as np
    
    # 指定默认字体
    pylab.rcParams['font.sans-serif'] = ['KaiTi']
    # 解决保存图像是负号'-'显示为方块的问题
    pylab.rcParams['axes.unicode_minus'] = False
    im = np.array(Image.open(r'D:\image_processing\image3\q1.jpg').convert('L'))
    freq = fp.fft2(im)
    (w, h) = freq.shape
    half_w, half_h = int(w/2), int(h/2)
    freq1 = np.copy(freq)
    freq2 = fftpack.fftshift(freq1)
    freq2_low = np.copy(freq2)
    freq2_low[half_w-10:half_w+11, half_h-10:half_h+11] = 0
    freq2 -= freq2_low
    im1 = fp.ifft2(fftpack.ifftshift(freq2)).real
    pylab.figure(figsize=(20,15))
    pylab.subplot(131)
    pylab.imshow((20*np.log10(0.1+freq2)).astype(int))
    pylab.title('LPF后的频谱', size=25)
    pylab.show()

    image.png



    以下代码演示了在不同的截止频率F下,LPF灰度图像上的应用.

    from PIL import Image
    import scipy.fftpack as fp
    from scipy import stats, fftpack
    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\image3\w2.jpg').convert('L'))
    freq = fp.fft2(im)
    (w, h) = freq.shape
    half_w, half_h = int(w/2), int(h/2)
    snrs_lp = []
    ubs = list(range(1, 25))
    pylab.figure(figsize=(12, 20))
    for u in ubs:
        freq1 = np.copy(freq)
        freq2 = fftpack.fftshift(freq1)
        freq2_low = np.copy(freq2)
        freq2_low[half_w-u:half_w+u+1,half_h-u:half_w+u+1] = 0
        freq2 -= freq2_low
        im1 = fp.ifft2(fftpack.ifftshift(freq2)).real
        snrs_lp.append(signaltonoise(im1))
        pylab.subplot(6, 4, u)
        pylab.imshow(im1, cmap='gray')
        pylab.axis('off')
        pylab.title('F = ' + str(u), size=20)
    pylab.subplots_adjust(wspace=0.1, hspace=0)
    pylab.show()

    image.png

    image.png

    可以看到,随着截止频率F的增大,LPF检测到的图像细节越丰富。



关键字