精华内容
下载资源
问答
  • 频率域滤波

    千次阅读 2019-11-11 19:38:34
    之前介绍了关于空间域滤波的一些操作,主要目的是减少噪声和平滑图像,空间域的处理方法都是...频率域滤波(即波数)为图像像元的灰度值随着位置变化的空间频率,以频谱表示信息分布特征,可以将一幅图像像元值在空间...

    之前介绍了关于空间域滤波的一些操作,主要目的是减少噪声和平滑图像,空间域的处理方法都是在已知噪声数据类型基础上进行噪声过滤的,而在某些噪声图像类别上上述方法就显得十分渺小了,频率域的引入就是来寻找噪声的类型,并进行空间过滤的。其中频率是单位时间变化期间,一个周期函数重复相同值序列的次数。频率域滤波(即波数)为图像像元的灰度值随着位置变化的空间频率,以频谱表示信息分布特征,可以将一幅图像像元值在空间上的变化分解为 具有不同振幅、空间频率和相位的简振函数的线性叠加,图像中各种空间频率成分的组成和分布称为空间频谱。这种对图像的空间频率特征进行分解、处理和分析称为频率域处理。频率中最典型的就是傅立叶变换,傅立叶变换能把图像从空间域变换到只包含不同频率信息的频率域,原图像上的灰度突变部位、图像结构复杂的区域、图像细节及干扰噪声等信息集中在高频区,而原图像上灰度变化平缓部位的信息集中在低频区(变换函数还有小波变换,将在后续进行分析和解析)。

    傅立叶变换

    傅立叶变换是将时间域转换为频率域的工具,傅立叶指出任何周期或函数都可以表示为不同频率的正弦和/或余弦之和的形式,每个正弦和/或余弦乘以不同的系数(称该和为傅立叶级数)。无论函数多么复杂,只要它是周期的,并且满足一定条件,都可以用这样的和表示。一维变换函数(一对,变换及反变化函数)为:

    其中,变换后的的函数F(x),f(t)为变换前周期函数,μ为周期函数分量,周期从0到1/2π

    利用一维变量扩展到二维变量,二维傅立叶变换公式为:

    e^j2π(xμ+yv)实际上是一个平面波(的共轭),它的传播方向与(u,v)相同,角频率是 根号下x^2+y^2。可以看做是平面波在x方向的角频率是u和y方向的角频率是v。二维傅立叶变换具有旋转不变性了(证明:如何证明二维离散傅里叶变换的旋转不变性?)。同时,变换完,中心是两个角频率都为0(直流分量),越靠近中心,两个方向的角频率越低,合成的角频率与到中心的距离成正比,所以中心是低频分量,外部是高频分量。幅度决定图像的强弱,相位觉得图像的像素。分析

    频率域滤波基本步骤

    preprocessing:图像预处理步骤,用0填充f(x,y),在二维循环卷积时为了保证周期性进行0填充(傅里叶变换滤波时,为什么需要对输入数据进行零填充?

    Fourier transform:图像傅立叶变换

    filter function:滤波器

    Inverse Fourier transform:反向傅立叶变换

    postprocessing:图像裁剪,输出g(x,y)

    伪流程:

    1.f(x,y)是由MxN构成,填充图像大小为P(P=2M)xQ(Q=2N),f_p(x,y)

    2.将f_p(x,y)利用-1^(x+y)将图片滤波器移动到中心,低频部分移到中间,高频部分移到四周,以便后面的计算。

    3.计算f_p(x,y)傅立叶变化得到F(u,v)

    4.生成一个实的、对阵的过滤器,中心在(P/2,Q/2),大小为PxQ

    5.计算H(u,v)F(u,v),并进行换位,得到初始时图像的傅立叶变换(DEF)的数据排列形式,即低频部分在四周,高频部分在中间(原点在左上角意味着低频部分在左上角,又因为的DEF是中心对称的,可得初始时图像的DEF应是低频在四周高频在中间)。

    6.计算反傅立叶函数,得到每个元素的幅度(实部和虚部平方和的平方根),并取左上角MxN区域,得到最终的输出图像g(x,y)。

    过滤器

    将一幅图像的模糊化(频率域平滑)可通过高频的衰减来达到,称为低通滤波器。

    将一幅图像的锐利化(边缘聚变与高频有关)可通过低频的衰减来达到,称为高通滤波器。

    a.低通滤波器

    (1)理想低通滤波器

    在以原点为圆心,以D_{0}为半径的园内,无衰减地通过所有频率,而在该圆外"切断"所有频率的二维低通滤波器,称为理想低通滤波器,它由下面函数确定:

    H(u,v)=\left\{\begin{matrix} 1 & D(u,v)<=D_{0} \\ 0 & D(u,v) > D_{0} \end{matrix}\right.

    其中,D(u,v)是频率域中点(u,v)与频率矩阵中心的距离,即:

    D(u,v)=[(u-P/2)^{2}+(v-Q/2)^{^{2}}]^{1/2}

    图4.40(a)中是H(u,v)的透视图,图b显示了以图像显示的滤波器,表明在半径D0的圆内,所有频率都无衰减的通过,而在此圆外的所有频率则是完全被衰减(滤除)。理想低通滤波器是关于原点径向对称的,这意味着该滤波器完全是一个径向横截面来定义,将图绕着H(u,v) 360度旋转则得到图a。

    (2)布特沃斯低通滤波器

    截止频率位于距原点D0处的n阶布特沃斯低通滤波器(BLPF)的传递函数定义如下:
    H(u,v)=\frac{1}{1+[D(u,v)/D_{0}]^{2n}}

    BLPF传递函数并没有在通过频率和滤除频率之间给出明显截止的尖锐的不连续性。图4.44(a)是透视图,表明在半径D0的圆内,频率是有规则的通过,而在此圆外的所有频率则是完全被衰减(滤除)。在空间域的一阶布特沃斯滤波器没有振铃现象,在二阶滤波器中,振铃现象通过很难察觉(通过设置不同的n值,二阶的BLPF是在有效的低通滤波和可接受的振铃特性之间好的折中),单更高阶数的滤波器中振铃现象会很明显。

    (3)高斯低通滤波器

    二维高斯低通滤波器(GLPF)形式如下:

    H(u,v)=e^{-D^{2}(u,v)/2D_{0}^{2}}

    其中D(u,v)是距离频率矩形中心的距离,D_{0}是关于中心扩展度的度量。图4.47显示了一个GLPF函数的透视图、图像显示和径向剖面图。GLPF的傅立叶反变换也是高斯的,从图a看得到的空间高斯滤波器没有振铃。

    b.高通滤波器

    通过衰减图像的傅立叶变换的高频成分可以平滑图像。因为边缘和其他灰度的急剧变化与高频分量有关,所以图像的锐化可在频率域通过高通滤波来实现,高通滤波会衰减傅立叶变换中的低频分量而不会扰乱高频信息。

    (1)理想高通滤波器

    一个二维理想高通滤波器(IHPF)定义为

    H(u,v)=\left\{\begin{matrix} 0 & D(u,v)<=D_{0} \\ 1 & D(u,v) > D_{0} \end{matrix}\right.

    其中D_{0}是截止频率,D(u,v)是频率域中点(u,v)与频率矩阵中心的距离。在这种滤波器下,IHPF把以半径为D_{0}的圆内的所有频率置0,而毫无衰减地通过圆外的所有频率。和理想低通过滤器完全相反。下图为理想高通滤波器的透视图、图像表示和剖面图。

    (2)布特沃斯高通滤波器

    截止频率D_{0}的n阶布特沃斯高通滤波器(BHPF)定义为:

    H(u,v)=\frac{1}{1+[D_{0}/D(u,v)]^{2n}}

    类似于低通滤波器,BHPF比IHPF更平滑一些。下图为布特沃斯高通滤波器的透视图、图像表示和剖面图。

    (3)高斯高通滤波器

    截止频率出在距频率矩形中心距离为D_{0}的高斯高通滤波器(GHPF)的传递函数:

    H(u,v)=e^{-2D_{0}^{2}/D^{2}(u,v)}

    类似高斯低通滤波器,高斯高通滤波器比理想高通滤波器、布特沃斯高通滤波器平滑效果更好一些。即使是对微小物体和细线条使用高斯滤波器滤波,结果也是比较清晰的。下图为高斯高通滤波器的透视图、图像表示和剖面图。

     

    https://www.cnblogs.com/laumians-notes/p/8544256.html

    https://blog.csdn.net/liuweiyuxiang/article/details/77040942

    https://blog.csdn.net/nanhuaibeian/article/details/90738701

    https://blog.csdn.net/chanxiaogang/article/details/45226373

    https://blog.csdn.net/eeeee123456/article/details/82950986

    展开全文
  • 目录频率域滤波基础频率域的其他特性频率域滤波基础知识频率域滤波步骤小结空间域和频率域滤波之间的对应关系 频率域滤波基础 频率域的其他特性 频率域中的滤波过程如下: 首先修改傅里叶变换以在到特定目的 然后...

    频率域滤波基础

    频率域的其他特性

    • 频率域中的滤波过程如下:
      • 首先修改傅里叶变换以在到特定目的
      • 然后计算IDFT,返回到空间域
    # 频率域中的其他特性
    img = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0429(a)(blown_ic).tif', -1)
    
    # FFT
    img_fft = np.fft.fft2(img.astype(np.float32))
    
    # 非中心化的频谱
    spectrum = spectrum_fft(img_fft)
    spectrum = np.uint8(normalize(spectrum) * 255)
    
    # 中心化
    fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心
    
    # 中心化后的频谱
    spectrum_fshift = spectrum_fft(fshift)
    spectrum_fshift_n = np.uint8(normalize(spectrum_fshift) * 255)
    
    # 对频谱做对数变换
    spectrum_log = np.log(1 + spectrum_fshift)
    
    plt.figure(figsize=(15, 8))
    plt.subplot(1, 2, 1), plt.imshow(img, cmap='gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(1, 2, 2), plt.imshow(spectrum_log, cmap='gray'), plt.title('FFT Shift log'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    频率域滤波基础知识

    已知大小为 P × Q P \times Q P×Q像素的一幅数字图像 f ( x , y ) f(x, y) f(x,y),则我们感兴趣取的基本滤波公式是:
    g ( x , y ) = R e a l { J − 1 [ H ( u , v ) F ( u , v ) ] } (4.104) g(x, y) = Real\{ \mathfrak{J}^{-1}[H(u, v) F(u, v)] \} \tag{4.104} g(x,y)=Real{J1[H(u,v)F(u,v)]}(4.104)

    J − 1 \mathfrak{J}^{-1} J1是IDFT, H ( u , v ) H(u, v) H(u,v)是滤波器传递函数(经常称为滤波器或滤波器函数), F ( u , v ) F(u, v) F(u,v)是输入图像 f ( x , y ) f(x, y) f(x,y)的DFT, g ( x , y ) g(x, y) g(x,y)是滤波后的图像。

    中以化:用 ( − 1 ) x + y (-1)^{x + y} (1)x+y乘以输入图像来完成的

    • 低通滤波器

      • 会衰减高频而通过低频的函数,会模糊图像
    • 高通滤波器

      • 将使图像的细节更清晰,但会降低图像的对比度
    def centralized_2d_loop(arr):
        """
        centralized 2d array $f(x, y) (-1)^{x + y}$
        """    
        #==================中心化=============================
        height, width = arr.shape[:2]
        dst = np.zeros([height, width])
        
        for h in range(height):
            for w in range(width):
                dst[h, w] = arr[h, w] * ((-1) ** (h + w))
        return dst
    
    def centralized_2d_index(arr):
        """
        centralized 2d array $f(x, y) (-1)^{x + y}$, about 35 times faster than loop
        """
        
        def get_1_minus1(img):
            """
            get 1 when image index is even, -1 while index is odd, same shape as input image, need this array to multiply with input image
            to get centralized image for DFT
            Parameter: img: input, here we only need img shape to create the array
            return such a [[1, -1, 1], [-1, 1, -1]] array, example is 3x3
            """
            M, N = img.shape
            m, n = img.shape
            # 这里判断一个m, n是否是偶数,如果是偶数的话,结果不正确,需要变换为奇数
            if m % 2 == 0 :
                m += 1
            if n % 2 ==0 :
                n += 1
            temp = np.arange(m * n).reshape(m, n)
            dst = np.piecewise(temp, [temp % 2 == 0, temp % 2 == 1], [1, -1])
            dst = dst[:M, :N]
            return dst
    
        #==================中心化=============================
        mask = get_1_minus1(arr)
        dst = arr * mask
    
        return dst
    
    def centralized_2d(arr):
        """
        centralized 2d array $f(x, y) (-1)^{x + y}$, about 4.5 times faster than index, 160 times faster than loop,
        """
        
        def get_1_minus1(img):
            """
            get 1 when image index is even, -1 while index is odd, same shape as input image, need this array to multiply with input image
            to get centralized image for DFT
            Parameter: img: input, here we only need img shape to create the array
            return such a [[1, -1, 1], [-1, 1, -1]] array, example is 3x3
            """
            dst = np.ones(img.shape)
            dst[1::2, ::2] = -1
            dst[::2, 1::2] = -1
            return dst
    
        #==================中心化=============================
        mask = get_1_minus1(arr)
        dst = arr * mask
    
        return dst
    
    a = np.arange(25).reshape(5, 5)
    a
    
    array([[ 0,  1,  2,  3,  4],
           [ 5,  6,  7,  8,  9],
           [10, 11, 12, 13, 14],
           [15, 16, 17, 18, 19],
           [20, 21, 22, 23, 24]])
    
    centralized_2d(a)
    
    array([[  0.,  -1.,   2.,  -3.,   4.],
           [ -5.,   6.,  -7.,   8.,  -9.],
           [ 10., -11.,  12., -13.,  14.],
           [-15.,  16., -17.,  18., -19.],
           [ 20., -21.,  22., -23.,  24.]])
    
    def idea_high_pass_filter(source, radius=5):
        M, N = source.shape[1], source.shape[0]
        
        u = np.arange(M)
        v = np.arange(N)
        
        u, v = np.meshgrid(u, v)
        
        D = np.sqrt((u - M//2)**2 + (v - N//2)**2)
        D0 = radius
        
        kernel = D.copy()
        
        kernel[D > D0] = 1
        kernel[D <= D0] = 0
        
        kernel_normal = (kernel - kernel.min()) / (kernel.max() - kernel.min())
        return kernel_normal
    
    # 频率域滤波基础知识
    img = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0429(a)(blown_ic).tif', -1)
    
    # FFT
    img_fft = np.fft.fft2(img.astype(np.float32))
    
    # 中心化
    fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心
    
    # 滤波器
    h = idea_high_pass_filter(img, radius=1)
    
    # 滤波
    res = fshift * h
    
    # 先反中以化,再反变换
    f1shift = np.fft.ifftshift(res)
    ifft = np.fft.ifft2(f1shift)
    
    # 取实数部分
    recon = np.real(ifft)
    # 小于0都被置为0
    # recon = np.where(recon > 0, recon, 0)
    recon = np.clip(recon, 0, recon.max())
    recon = np.uint8(normalize(recon) * 255)
    
    plt.figure(figsize=(15, 8))
    plt.subplot(1, 2, 1), plt.imshow(img, cmap='gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(1, 2, 2), plt.imshow(recon, cmap='gray'), plt.title('FFT Shift log'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    def gauss_high_pass_filter(source, center, radius=5):
        """
        create gaussian high pass filter 
        param: source: input, source image
        param: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is 
                       center of the source image
        param: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then all
                       value is 1
        return a [0, 1] value filter
        """
        M, N = source.shape[1], source.shape[0]
        
        u = np.arange(M)
        v = np.arange(N)
        
        u, v = np.meshgrid(u, v)
        
        D = np.sqrt((u - center[1]//2)**2 + (v - center[0]//2)**2)
        D0 = radius
        
        if D0 == 0:
            kernel = np.ones(source.shape[:2], dtype=np.float)
        else:
            kernel = 1 - np.exp(- (D**2)/(D0**2))   
            kernel = (kernel - kernel.min()) / (kernel.max() - kernel.min())
        return kernel
    
    def plot_3d(ax, x, y, z):
        ax.plot_surface(x, y, z, antialiased=True, shade=True)
        ax.view_init(30, 60), ax.grid(b=False), ax.set_xticks([]), ax.set_yticks([]), ax.set_zticks([])
        
    def imshow_img(img, ax, cmap='gray'):
        ax.imshow(img, cmap=cmap), ax.set_xticks([]), ax.set_yticks([])
    
    def frequency_filter(fshift, h):
        """
        Frequency domain filtering using FFT
        param: fshift: input, FFT shift to center
        param: h: input, filter, value range [0, 1]
        return a uint8 [0, 255] reconstruct image
        """
        # 滤波
        res = fshift * h
    
        # 先反中以化,再反变换
        f1shift = np.fft.ifftshift(res)
        ifft = np.fft.ifft2(f1shift)
    
        # 取实数部分
        recon = np.real(ifft)
        # 小于0都被置为0
        # recon = np.where(recon > 0, recon, 0)
        recon = np.clip(recon, 0, recon.max())
        recon = np.uint8(normalize(recon) * 255)
        
        return recon
    
    # 高斯核滤波器与对应的结果
    from mpl_toolkits.mplot3d import Axes3D
    import numpy as np
    from matplotlib import pyplot as plt
    from matplotlib import cm
    
    # 显示多个3D高斯核函数
    k = 3
    sigma = 0.8
    X = np.linspace(-k, k, 200)
    Y = np.linspace(-k, k, 200)
    x, y = np.meshgrid(X, Y)
    z = np.exp(- ((x)**2 + (y)**2)/ (2 * sigma**2))
    
    fig = plt.figure(figsize=(15, 10))
    ax_1 = fig.add_subplot(2, 3, 1, projection='3d')
    plot_3d(ax_1, x, y, z)
    
    ax_2 = fig.add_subplot(2, 3, 2, projection='3d')
    plot_3d(ax_2, x, y, 1 - z)
    
    # 这里只显示不明显,
    a = 1
    x = a + x
    y = a + y
    ax_3 = fig.add_subplot(2, 3, 3, projection='3d')
    plot_3d(ax_3, x, y, 1 - z)
    
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)
    
    # FFT
    img_fft = np.fft.fft2(img_ori.astype(np.float32))
    
    # 中心化
    fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心
    
    # 滤波器
    radius = 20
    center = img_ori.shape[:2]
    h_high_1 = gauss_high_pass_filter(img_ori, center, radius=radius)
    h_low = 1 - h_high_1
    center = (1026, 872)
    h_high_2 = gauss_high_pass_filter(img_ori, center, radius=radius)
    
    img_low = frequency_filter(fshift, h_low)
    img_high_1 = frequency_filter(fshift, h_high_1)
    img_high_2 = frequency_filter(fshift, h_high_2)
    
    ax_4 = fig.add_subplot(2, 3, 4)
    imshow_img(img_low, ax_4)
    
    ax_5 = fig.add_subplot(2, 3, 5)
    imshow_img(img_high_1, ax_5)
    
    ax_6 = fig.add_subplot(2, 3, 6)
    imshow_img(img_high_2, ax_6)
    
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    下面展示的是函数未填充时,会产生交叠错误;而函数填充0后,产生我们期望的结果,但是此时的滤波器应该是原来的2倍大小,如果是还是原来的大小的话,那么图像会比之前更模糊的。

    这里采用的方法是对图像进行零填充,再直接在频率域中创建期望的滤波器传递函数,该函数的到底会给你个填充零后的图像的大小相同。这个明显降低交叠错误,只是会出现振铃效应

    # 未填充与填充0的结果
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0432(a)(square_original).tif', -1)
    
    #================未填充===================
    # FFT
    img_fft = np.fft.fft2(img_ori.astype(np.float32))
    
    # 中心化
    fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心
    
    # 滤波器
    radius = 50
    center = img_ori.shape[:2]
    h_high_1 = gauss_high_pass_filter(img_ori, center, radius=radius)
    h_low = 1 - h_high_1
    
    img_low = frequency_filter(fshift, h_low)
    
    fig = plt.figure(figsize=(15, 5))
    ax_1 = fig.add_subplot(1, 3, 1)
    imshow_img(img_ori, ax_1)
    
    ax_2 = fig.add_subplot(1, 3, 2)
    imshow_img(img_low, ax_2)
    
    #================填充0====================
    fp = np.zeros([img_ori.shape[0]*2, img_ori.shape[1]*2])
    fp[:img_ori.shape[0], :img_ori.shape[1]] = img_ori
    
    img_fft = np.fft.fft2(fp.astype(np.float32))
    fshift = np.fft.fftshift(img_fft)   
    
    hp = 1 - gauss_high_pass_filter(fp, fp.shape, radius=100)
    
    img_low = frequency_filter(fshift, hp)
    img_dst = img_low[:img_ori.shape[0], :img_ori.shape[1]]
    
    ax_3 = fig.add_subplot(1, 3, 3)
    imshow_img(img_dst, ax_3)
    
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    另一个合理的做法是:构建一个与未填充图像大小相同的函数,计算该函数的IDFT得到相应的空间表示,在空间域中对这一表示填充,然后计算其DFT返回到频率域。这也会出现振铃效应

    def centralized(arr):
        """
        centralized 1d array $f(x) (-1)^{n}$
        """
        u = arr.shape[0]
        dst = np.zeros(arr.shape)
        
        # 中心化
        for n in range(u):
            dst[n] = arr[n] * ((-1) ** n)
        return dst
    
    # 频率域反变换到空间域,再填充,再DFT,会出现振铃效应
    h = np.zeros([256])
    h[:9] = 1
    h_c = np.fft.fftshift(h)
    
    fig = plt.figure(figsize=(15, 10))
    grid = plt.GridSpec(2, 3, wspace=0.2, hspace=0.1)
    ax_1 = fig.add_subplot(grid[0, 0])
    ax_1.plot(h_c), ax_1.set_xticks([0, 128, 255]), ax_1.set_yticks(np.arange(-0.2, 1.4, 0.2)), ax_1.set_xlim([0, 255])
    
    # 用公式法的中以化才能得到理想的结果
    h_pad = np.pad(h, (0, 256), mode='constant', constant_values=0)
    ifshift = centralized_2d(h_pad.reshape(1, -1)).flatten() # 也可以用centralized_2d这个函数
    ifft = np.fft.ifft(ifshift)
    ifft = ifft.real * 2
    ifft[:128] = 0
    ifft[384:] = 0
    ax_2 = fig.add_subplot(grid[0, 1:3])
    ax_2.plot(ifft), ax_2.set_xticks([0, 128, 256, 384, 511]), ax_2.set_yticks(np.arange(-0.01, 0.05, 0.01)), ax_2.set_xlim([0, 511])
    
    f = ifft[128:384]
    ax_3 = fig.add_subplot(grid[1, 0])
    ax_3.plot(f), ax_3.set_xticks([0, 128, 255]), ax_3.set_yticks(np.arange(-0.01, 0.05, 0.01)), ax_3.set_xlim([0, 255])
    
    # 已经中心化,截取部分构成原函数,但出现的效果不是太好,但可以看出振铃效应
    d = 120
    f_ = np.zeros([h.shape[0]])
    f_[:256-d] = f[d:]
    f_pad = np.pad(f_, (0, 256), mode='constant', constant_values=0)
    f_centralized = centralized(f_pad)
    fft = np.fft.fft(f_centralized)
    ax_4 = fig.add_subplot(grid[1, 1:3])
    ax_4.plot(fft.real), ax_4.set_xticks([0, 128, 256, 384, 511]), ax_4.set_yticks(np.arange(-0.2, 1.4, 0.2)), ax_4.set_xlim([0, 511])
    
    plt.show()
    

    在这里插入图片描述

    零相移滤波器

    • 相角计算为复数的实部与虚部之比的反正切。因为 H ( u , v ) H(u, v) H(u,v)同是乘以 R , I R,I R,I,所以当这个比率形成后,它会被抵消。对实部和虚部的影响相同且对相角无影响的滤波器。
    # 相角对反变换的影响
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)
    
    #================未填充===================
    # FFT
    img_fft = np.fft.fft2(img_ori.astype(np.float32))
    
    # 中心化
    fshift = np.fft.fftshift(img_fft)            
    
    # 频谱与相角
    spectrum = spectrum_fft(fshift)
    phase = phase_fft(fshift)
    
    # 相角乘以-1
    phase_1 = phase * -1
    new_fshift = spectrum * np.exp(phase_1 * -1j)
    new_fshift = np.fft.ifftshift(new_fshift)
    recon_1 = np.fft.fft2(new_fshift)
    
    # 相角乘以0.5,需要乘以一个倍数,以便可视化
    phase_2 = phase * 0.25
    new_fshift = spectrum * np.exp(phase_2 * -1j)
    new_fshift = np.fft.ifftshift(new_fshift)
    recon_2 = np.fft.fft2(new_fshift)
    recon_2 = np.uint8(normalize(recon_2.real) * 4000)
    
    fig = plt.figure(figsize=(15, 5))
    ax_1 = fig.add_subplot(1, 3, 1)
    imshow_img(img_ori, ax_1)
    
    ax_2 = fig.add_subplot(1, 3, 2)
    imshow_img(recon_1.real, ax_2)
    
    ax_3 = fig.add_subplot(1, 3, 3)
    ax_3.imshow(recon_2, cmap='gray')
    
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    频率域滤波步骤小结

    1. 已知一幅大小为 M × N M\times N M×N的输入图像 f ( x , y ) f(x, y) f(x,y),分别得到填充零后的尺寸 P P P Q Q Q,即 P = 2 M P=2M P=2M Q = 2 N Q=2N Q=2N
    2. 使用零填充、镜像填充或复制填充,形成大小为 P × Q P\times Q P×Q的填充后的图像 f P ( x , y ) f_P(x, y) fP(x,y)
    3. f P ( x , y ) f_P(x, y) fP(x,y)乘以 ( − 1 ) x + y (-1)^{x+y} (1)x+y,使用傅里叶变换位于 P × Q P\times Q P×Q大小的频率矩形的中心。
    4. 计算步骤3得到的图像的DFT,即 F ( u , v ) F(u, v) F(u,v)
    5. 构建一个实对称滤波器传递函数 H ( u , v ) H(u, v) H(u,v),其大小为 P × Q P\times Q P×Q,中心中 ( P / 2 , Q / 2 ) (P/2, Q/2) (P/2,Q/2)处。
    6. 采用对应像素相乘得到 G ( u , v ) = H ( u , v ) F ( u , v ) G(u, v) = H(u, v)F(u, v) G(u,v)=H(u,v)F(u,v)
    7. 计算 G ( u , v ) G(u, v) G(u,v)的IDFT得到滤波后的图像(大小为 P × Q P\times Q P×Q), g P ( x , y ) = { Real [ J − 1 [ G ( u , v ) ] ] } ( − 1 ) x + y g_P(x, y) = \Big\{\text{Real} \big[\mathfrak{J}^{-1}[G(u, v)] \big]\Big\}(-1)^{x+y} gP(x,y)={Real[J1[G(u,v)]]}(1)x+y
    8. 最后,从 g P ( x , y ) g_P(x, y) gP(x,y)的左上象限提取一个大小为 M × N M\times N M×N的区域,得到与输入图像大小相同的滤波后的结果 g ( x , y ) g(x, y) g(x,y)
    def pad_image(img, mode='constant'):
        """
        pad image into PxQ shape, orginal is in the top left corner.
        param: img: input image
        param: mode: input, str, is numpy pad parameter, default is 'constant'. for more detail please refere to Numpy pad
        return PxQ shape image padded with zeros or other values
        """
        dst = np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), mode=mode)
        return dst    
    
    # 频率域滤波过程
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)
    M, N = img_ori.shape[:2]
    fig = plt.figure(figsize=(15, 15))
    
    # 这里对原图像进行pad,以得到更好的显示
    img_ori_show = np.ones([img_ori.shape[0]*2, img_ori.shape[1]*2]) * 255
    img_ori_show[:img_ori.shape[0], img_ori.shape[1]:] = img_ori
    ax_1 = fig.add_subplot(3, 3, 1)
    imshow_img(img_ori_show, ax_1)
    
    fp = pad_image(img_ori, mode='constant')
    ax_2 = fig.add_subplot(3, 3, 2)
    imshow_img(fp, ax_2)
    
    fp_cen = centralized_2d(fp)
    fp_cen_show = np.clip(fp_cen, 0, 255) # 负数的都截断为0
    ax_3 = fig.add_subplot(3, 3, 3)
    ax_3.imshow(fp_cen_show, cmap='gray'), ax_3.set_xticks([]), ax_3.set_yticks([])
    
    fft = np.fft.fft2(fp_cen)
    spectrum = spectrum_fft(fft)
    spectrum = np.log(1 + spectrum)
    ax_4 = fig.add_subplot(3, 3, 4)
    imshow_img(spectrum, ax_4)
    
    H = 1 - gauss_high_pass_filter(fp, center=fp.shape, radius=70)
    ax_5 = fig.add_subplot(3, 3, 5)
    imshow_img(H, ax_5)
    
    HF = fft * H
    HF_spectrum = np.log(1 + spectrum_fft(HF))
    ax_6 = fig.add_subplot(3, 3, 6)
    imshow_img(HF_spectrum, ax_6)
    
    ifft = np.fft.ifft2(HF)
    gp = centralized_2d(ifft.real)
    ax_7 = fig.add_subplot(3, 3, 7)
    imshow_img(gp, ax_7)
    
    # 最终结果有黑边
    g = gp[:M, :N]
    ax_8 = fig.add_subplot(3, 3, 8)
    imshow_img(g, ax_8)
    
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    def frequency_filter(fft, h):
        """
        Frequency domain filtering using FFT
        param: fshift: input, FFT shift to center
        param: h: input, filter, value range [0, 1]
        return a uint8 [0, 255] reconstruct image
        """
        # 滤波
        res = fft * h
    
        # 先反中以化,再反变换
        ifft = np.fft.ifft2(res)
    
        # 取实数部分
        recon = centralized_2d(ifft.real)
        # 小于0都被置为0
        recon = np.clip(recon, 0, recon.max())
        recon = np.uint8(normalize(recon) * 255)
        
        return recon
    
    # 频率域滤波
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0431(d)(blown_ic_crop).tif', -1)
    M, N = img_ori.shape[:2]
    
    # 填充
    fp = pad_image(img_ori, mode='reflect')
    # 中心化
    fp_cen = centralized_2d(fp)
    # 正变换
    fft = np.fft.fft2(fp_cen)
    # 滤波器
    H = 1 - gauss_high_pass_filter(fp, center=fp.shape, radius=100)
    # 滤波
    HF = fft * H
    # 反变换
    ifft = np.fft.ifft2(HF)
    # 去中心化
    gp = centralized_2d(ifft.real)
    # 还回来与原图像的大小
    g = gp[:M, :N]
    fig = plt.figure(figsize=(15, 15))
    ax_8 = fig.add_subplot(3, 3, 8)
    imshow_img(g, ax_8)
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    空间域和频率域滤波之间的对应关系

    空间域中的滤波器核与频率域傅里叶变换对:
    h ( x , y ) ⇔ H ( u , v ) (4.106) h(x, y) \Leftrightarrow H(u, v) \tag{4.106} h(x,y)H(u,v)(4.106)

    h ( x , y ) h(x, y) h(x,y)是空间核。这个核 可以由一个频率域滤波器对一个冲激的响应得到,因此有时我们称 h ( x , y ) h(x, y) h(x,y) H ( u , v ) H(u, v) H(u,v)的冲激响应。

    一维频率域高斯传递函数
    H ( u ) = A e − u 2 / 2 σ 2 (4.107) H(u) = A e^{-u^2/2\sigma^2} \tag{4.107} H(u)=Aeu2/2σ2(4.107)
    空间域的核 由 H ( u ) H(u) H(u)的傅里叶反变换得到
    h ( x ) = 2 π σ A e − 2 π 2 σ 2 x 2 (4.108) h(x) = \sqrt{2\pi}\sigma A e^{-2\pi^2\sigma^2 x^2} \tag{4.108} h(x)=2π σAe2π2σ2x2(4.108)

    这是两个很重要的公式

    1. 它们是一个傅里叶变换对,两者都是高斯函数和实函数
    2. 两个函数的表现相反
    • 高斯差

      • 用一个常量减去低通函数,可以构建一个高通滤波器。
      • 因为使用所谓的高斯差(涉及两个低通函数)就可更好地控制滤波函数的形状
    • 频率域 A ≥ B , σ 1 > σ 2 A \ge B, \sigma_1 > \sigma_2 AB,σ1>σ2
      H ( u ) = A e − u 2 / 2 σ 1 2 − B e − u 2 / 2 σ 2 2 (4.109) H(u) = A e^{-u^2/2\sigma_1^2} - B e^{-u^2/2\sigma_2^2}\tag{4.109} H(u)=Aeu2/2σ12Beu2/2σ22(4.109)

    • 空间域
      h ( x ) = 2 π σ 1 A e − 2 π 2 σ 1 2 x 2 − 2 π σ 2 B e − 2 π 2 σ 2 2 x 2 (4.110) h(x) = \sqrt{2\pi}\sigma_1 A e^{-2\pi^2\sigma_1^2 x^2} - \sqrt{2\pi}\sigma_2 B e^{-2\pi^2\sigma_2^2 x^2} \tag{4.110} h(x)=2π σ1Ae2π2σ12x22π σ2Be2π2σ22x2(4.110)

    # 空间域与频率域的滤波器对应关系
    u = np.linspace(-3, 3, 200)
    
    A = 1
    sigma = 0.2
    Hu = A * np.exp(-np.power(u, 2) / (2 * np.power(sigma, 2)))
    
    x = u
    hx = np.sqrt(2 * np.pi) * sigma * A * np.exp(-2 * np.pi**2 * sigma**2 * x**2)
    
    fig = plt.figure(figsize=(10, 10))
    ax_1 = setup_axes(fig, 221)
    ax_1.plot(u, Hu), ax_1.set_title('H(u)', loc='center', y=1.05), ax_1.set_ylim([0, 1]), ax_1.set_yticks([])
    
    B = 1
    sigma_2 = 4
    Hu_2 = B * np.exp(-np.power(u, 2) / (2 * np.power(sigma_2, 2)))
    Hu_diff = Hu_2 - Hu
    
    ax_2 = setup_axes(fig, 222)
    ax_2.plot(x, Hu_diff), ax_2.set_title('H(u)', loc='center', y=1.05), ax_2.set_ylim([0, 1.1]), ax_2.set_yticks([])
    
    ax_3 = setup_axes(fig, 223)
    ax_3.plot(x, hx), ax_3.set_title('h(x)', loc='center', y=1.05), ax_3.set_ylim([0, 0.8]), ax_3.set_yticks([])
    
    hx_2 = np.sqrt(2 * np.pi) * sigma_2 * B * np.exp(-2 * np.pi**2 * sigma_2**2 * x**2)
    hx_diff = hx_2 - hx
    
    ax_4 = setup_axes(fig, 224)
    ax_4.plot(x, hx_diff), ax_4.set_title('h(x)', loc='center', y=1.05), ax_4.set_ylim([-1, 10]), ax_4.set_yticks([])
    
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    # 空间核得到频率域
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0438(a)(bld_600by600).tif', -1)
    
    fp = pad_zero(img_ori)
    fp_cen = centralized_2d(fp)
    fft = np.fft.fft2(fp_cen)
    spectrum = spectrum_fft(fft)
    spectrum = np.log(1 + spectrum)
    
    plt.figure(figsize=(18, 15))
    plt.subplot(1, 2, 1), plt.imshow(img, cmap='gray'), plt.title('No Log')
    plt.subplot(1, 2, 2), plt.imshow(spectrum, cmap='gray'), plt.title('FFT')
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述
    注:频率域的Sobel滤波器还没有解决

    # 频域
    center = (fp.shape[0], fp.shape[1] + 300)
    GHPF = gauss_high_pass_filter(fp, center, radius=200)
    GHPF = (1 - np.where(GHPF > 1e-6, GHPF, 1)) * -2
    
    center = (fp.shape[0], fp.shape[1] - 300)
    GHPF_1 = gauss_high_pass_filter(fp, center, radius=200)
    GHPF_1 = (1 - np.where(GHPF_1 > 1e-6, GHPF_1, 1)) * 2
    
    sobel_f = GHPF + GHPF_1
    
    plt.figure(figsize=(15, 5))
    
    plt.subplot(1, 3, 1), plt.imshow(GHPF, 'gray')
    plt.subplot(1, 3, 2), plt.imshow(GHPF_1, 'gray')
    plt.subplot(1, 3, 3), plt.imshow(sobel_f, 'gray')
    
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    # 频域
    fig = plt.figure(figsize=(15, 5))
    HF = fft * sobel_f
    HF_spectrum = np.log(1 + spectrum_fft(HF))
    ax_6 = fig.add_subplot(1, 3, 1)
    imshow_img(HF_spectrum, ax_6)
    
    ifft = np.fft.ifft2(HF)
    gp = centralized_2d(ifft.real)
    ax_7 = fig.add_subplot(1, 3, 2)
    imshow_img(gp, ax_7)
    
    # 最终结果有黑边
    g = gp[:M, :N]
    ax_8 = fig.add_subplot(1, 3, 3)
    imshow_img(g, ax_8)
    
    # plt.figure(figsize=(15, 5))
    
    # plt.subplot(1, 3, 1), plt.imshow(GHPF, 'gray')
    # plt.subplot(1, 3, 2), plt.imshow(GHPF_1, 'gray')
    # plt.subplot(1, 3, 3), plt.imshow(sobel_f, 'gray')
    
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    空间域的Sobel滤波器

    def kernel_seperate(kernel):
        """
        这里的分离核跟之前不太一样
        get kernel seperate vector if separable
        param: kernel: input kerel of the filter
        return two separable vector c, r
        """
        rank = np.linalg.matrix_rank(kernel)
    
        if rank == 1:
    
        #-----------------Numpy------------------
        # here for sobel kernel seperate
            c = kernel[:, -1]
            r = kernel[-1, :]
    
        else:
            print('Not separable')
    
    #     c = np.array(c).reshape(-1, 1)
    #     r = np.array(r).reshape(-1, 1)
        
        return c, r
    
    def separate_kernel_conv2D(img, kernel):
        """
        separable kernel convolution 2D, if 2D kernel not separable, then canot use this fuction. in the fuction. first need to
        seperate the 2D kernel into two 1D vectors.
        param: img: input image you want to implement 2d convolution with 2d kernel
        param: kernel: 2D separable kernel, such as Gaussian kernel, Box Kernel
        return image after 2D convolution
        """
        m, n = kernel.shape[:2]
        pad_h = int((m - 1) / 2)
        pad_w = int((n - 1) / 2)
        w_1, w_2 = kernel_seperate(kernel)
    #     w_1 = np.array([-1, 0, 1])
    #     w_2 = np.array([1, 2, 1])
        convovle = np.vectorize(np.convolve, signature='(x),(y)->(k)')
        r_2 = convovle(img, w_2)
        r_r = np.rot90(r_2)
        r_1 = convovle(r_r, w_1)
        r_1 = r_1.T
        r_1 = r_1[:, ::-1]
        img_dst = img.copy().astype(np.float)
        img_dst = r_1[pad_h:-pad_h, pad_w:-pad_w]
        return img_dst
    
    # Sobel梯度增强边缘
    img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH04/Fig0438(a)(bld_600by600).tif", 0)
    
    sobel_x = np.zeros([3, 3], np.int)
    sobel_x[0, :] = np.array([-1, -2, -1])
    sobel_x[2, :] = np.array([1, 2, 1])
    
    sobel_y = np.zeros([3, 3], np.int)
    sobel_y[:, 0] = np.array([-1, -2, -1])
    sobel_y[:, 2] = np.array([1, 2, 1])
    
    gx = separate_kernel_conv2D(img_ori, kernel=sobel_x)
    gy = separate_kernel_conv2D(img_ori, kernel=sobel_y)
    
    # gx = cv2.filter2D(img_ori, ddepth=-1, kernel=sobel_x)
    # gy = cv2.filter2D(img_ori, ddepth=-1, kernel=sobel_y)
    
    # 先对gx gy做二值化处理再应用下面的公式
    img_sobel = np.sqrt(gx**2 + gy**2)   
    
    # 二值化后,平方根的效果与绝对值很接近
    img_sobel = abs(gx) + abs(gy)
    img_sobel = np.uint8(normalize(img_sobel) * 255)
    
    plt.figure(figsize=(15, 12))
    plt.subplot(1, 2, 1), plt.imshow(img_ori,   'gray', vmax=255), plt.title("Original"), plt.xticks([]), plt.yticks([])
    plt.subplot(1, 2, 2), plt.imshow(img_sobel, 'gray', vmax=255), plt.title("Sobel"), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    展开全文
  • 1.频率域滤波基础 1.基础知识 变化最慢的频率分量(u=v=0)与图像的平均灰度成正比;当我们远离变换的原点是,低频对应于图像中变化缓慢的灰度分量。变换幅度(谱)和相角。 频率域滤波:修改一幅图像的傅里叶变换、...

    版权声明:以下图片截取自《数字图像处理》冈萨雷斯 一书中。

    1.频率域滤波基础

    1.基础知识
    变化最慢的频率分量(u=v=0)与图像的平均灰度成正比;当我们远离变换的原点是,低频对应于图像中变化缓慢的灰度分量。变换幅度(谱)和相角
    频率域滤波:修改一幅图像的傅里叶变换、然后计算其反变换得到处理后的结果;其基本滤波公式如下(中心化):
    在这里插入图片描述
    变换中的低频:与图像中缓慢变化的灰度分量有关(室内的墙面和室外少云的天空等);
    变换中的高频:由灰度的尖锐过渡造成(边缘和噪声等)。
    低通滤波器(衰减高频而通过低频):模糊一幅图像;
    高通滤波器(衰减低频而通过高频):增强尖锐的细节,但将导致图像对比度的降低。
    DFT是复数阵列,可以将它表示为实部和虚部:
    在这里插入图片描述
    零相移滤波器:等同地影响实部和虚部而不影响相位的滤波器。
    频率域滤波步骤:
    在这里插入图片描述
    在这里插入图片描述
    2.空间滤波和频率域滤波
    空间域和频率域滤波间的纽带:卷积定理。
    给定一个空间滤波器,可以用该空间滤波器的傅里叶正变换得到其平吕与表示(傅里叶变换对):
    在这里插入图片描述
    高斯滤波器:
    一个高斯函数的正、反傅里叶变换都是实高斯函数。
    一维频率域高斯滤波器:
    在这里插入图片描述
    使用一个全部带正系数的模板就可以在空间域中实现低通滤波;
    滤波器宽度间的互易关系:频率域滤波器越窄,其衰减的低频越多,引起的模糊越大;在空间域,必须使用较大的模板来增加模糊。
    在这里插入图片描述

    2.使用频率域滤波器平滑图像

    低通滤波器:一幅图像中的边缘和其他尖锐的灰度转变(如,噪声)对其DFT的高频内容有贡献;因此,在频率域平滑(模糊)可通过对高频的衰减来达到。
    1.理想低通滤波器(ILPF)
    其公式如下:
    在这里插入图片描述
    其图示如下:
    在这里插入图片描述
    功率表示:
    在这里插入图片描述
    2.布特沃斯低通滤波器(BLPF)
    n阶布特沃斯低通滤波器:
    当阶数较高时,BLPF接近于理想滤波器;对于较低的阶数值,BLPF更像高斯滤波器。
    在这里插入图片描述
    在空间域的一阶布特沃斯滤波没有振铃现象;在二阶滤波器中,振铃现象通常很难察觉,但更高阶的滤波器中振铃现象会很明显。
    二阶的BLPF是在有效的低通滤波和可接受的振铃特性之间的好的折中。
    在这里插入图片描述
    3.高斯低通滤波器(GLPF)
    二维高斯滤波器:
    在这里插入图片描述
    GLPF的傅里叶反变换也是高斯的,这意味着通过计算式(4.8-6)和式(4.8-7)的IDFT得到的空间高斯滤波器将没有振铃。
    二维GLPF的图示如下:
    在这里插入图片描述

    4.小结
    表4.4:三类低通滤波器
    在这里插入图片描述

    3.使用频率域滤波器锐化图像

    边缘和其他灰度的急剧变化与高频分量有关,所以图像的锐化可在频率域通过高通滤波来实现;高通滤波会衰减DFT中的低频分量而不会扰乱高频信息。
    在这里插入图片描述
    1.理想高通滤波器(IHPF)、布特沃斯高通滤波器(BHPF)、高斯高通滤波器(GFHPF)
    IHPF:
    在这里插入图片描述
    BHPF:
    在这里插入图片描述
    GHPF:
    在这里插入图片描述
    图示:在这里插入图片描述
    2.频率域的拉普拉斯算子
    拉普拉斯算子在频率域的滤波实现:
    在这里插入图片描述
    其中,c=-1,因为H(u,v)时负的。
    在这里插入图片描述
    3.同态滤波
    关键:照射分量和反射分量的分离;
    在这里插入图片描述
    照射分量——慢的变化空间——DFT的低频;
    反射分量——突变、边缘信息——DFT的高频;
    在这里插入图片描述
    同态滤波用于控制照射分量和反射分量:
    在这里插入图片描述

    4.选择性滤波

    1.带阻滤波器和带通滤波器
    兴趣:处理指定的频段或频率矩形的小区域。
    D0:带宽的径向中心;
    W:带宽。
    带阻滤波器:
    在这里插入图片描述
    带通滤波器:
    在这里插入图片描述

    2.陷波滤波器
    陷波滤波器:拒绝(或通过)事先定义的关于频率矩形中心的一个领域的频率(更有用的选择性滤波器)。
    陷波带阻滤波器:用中心已被平移到陷波滤波器中心的高通滤波器的乘积来构造。一般形式如下:
    在这里插入图片描述
    陷波带通滤波器:
    在这里插入图片描述

    5.实现

    1.二维DFT的可分性
    f(x,y)的二维DFT可通过计算f(x,y)的每一行的一维变换,然后,沿着计算结果的每一列计算一维变换来得到。
    在这里插入图片描述
    2.用DFT算法来计算IDFT
    在这里插入图片描述

    展开全文
  • 使用opencv进行频率域滤波的实战,包括高通,低通,高斯,巴特沃斯等。
  • 频率域滤波示例

    2013-12-12 15:41:12
    频率域滤波,帮助理解信号处理课程中滤波问题,为matlab脚本程序
  • 有实现二维fft变换ifft变换方法以及各种频率域滤波以及实验用的图片和结果图
  • dft2d函数为对灰度图进行离散傅里叶变换和反变换 filter2d_freq为对灰度图进行频率域滤波 修改Runner函数中的图片路径然后运行即可
  • 数字图像处理——第4章 频率域滤波 文章目录数字图像处理——第4章 频率域滤波频率域1.傅里叶级数原理1.1.一维傅里叶变换1.2.二维傅里叶变换2.python×傅里叶级数2.1.傅里叶变换后的频谱图3.频率域滤波3.1.低频与...

    数字图像处理——第4章 频率域滤波

    频率域

    上一章学的灰度变换和空间滤波,主要目的是减少噪声和平滑图像,同时也是在空间域进行的图像增强操作。而这章主要是频率域滤波,滤波的概念在上一章节有涉及到,所以先了解下频率域。频率域是指从函数的频率角度出发分析函数,和频率域相对的是时间域。简单说就是如果从时间域分析信号时,时间是横坐标,振幅是纵坐标。而在频率域分析的时候则是频率是横坐标,振幅是纵坐标。那么下一个问题就是为什么要在频率域中进行图像处理? 1). 可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务,在频率域中变得非常普通;2). 滤波在频率域更为直观,它可以解释空间域滤波的某些性质; 3).可以在频率域指定滤波器,做反变换,然后在空间域使用结果滤波器作为空间域滤波器的指导。

    频率域滤波是对图像进行傅里叶变换,将图像由图像空间转换到频域空间,然后在频率域中对图像的频谱作分析处理,以改变图像的频率特征。

    1.傅里叶级数原理

    傅里叶是18世纪法国的一位伟大的数学家。他指出任何周期函数都可以表示为不同频率的正弦和或者余弦和的形式,每个正弦或者余弦乘以不同的系数就是傅里叶级数。无论函数有多复杂,只要它是周期性的,并且满足一定的数学条件,就一定可以用这样的正弦和或者余弦和的形式来表示。翻开考完研没扔的高数,回看到傅里叶级数,公式如下:
    f ( x ) = a 0 2 + ∑ n = 1 ∞ ( a n cos ⁡ n π l x + b n sin ⁡ n π l x ) f(x)=\frac{a_{0}}{2}+\sum_{n=1}^{\infty}\left(a_{n} \cos \frac{n \pi}{l} x+b_{n} \sin \frac{n \pi}{l} x\right) f(x)=2a0+n=1(ancoslnπx+bnsinlnπx)
    其中:

    a n = 1 l ∫ − l l f ( x ) cos ⁡ n π l x   d x , n = 0 , 1 , 2 , ⋯   , 其 中 a 0 为 n = 0 时 a_{n}=\frac{1}{l} \int_{-l}^{l} f(x) \cos \frac{n \pi}{l} x \mathrm{~d} x, n=0,1,2, \cdots,其中a_{0}为n=0时 an=l1llf(x)coslnπx dx,n=0,1,2,,a0n=0

    b n = 1 l ∫ − l l f ( x ) sin ⁡ n π l x   d x , n = 0 , 1 , 2 , ⋯   , b_{n}=\frac{1}{l} \int_{-l}^{l} f(x) \sin \frac{n \pi}{l} x \mathrm{~d} x, n=0,1,2, \cdots, bn=l1llf(x)sinlnπx dx,n=0,1,2,,

    除此之外,若展开式中只有正弦函数,则称其为正弦函数;相反,若展开式中只有余弦函数,则称其为余弦级数。

    1.1.一维傅里叶变换

    傅里叶变换:

    F ( μ ) = ∫ − ∞ + ∞ f ( t ) e − j 2 π μ t d t F(\mu)=\int_{-\infty}^{+\infty} f(t) e^{-j 2 \pi \mu t} d t F(μ)=+f(t)ej2πμtdt
    傅里叶反变换:
    f ( t ) = ∫ − ∞ + ∞ F ( μ ) e j 2 π μ t d μ f(t)=\int_{-\infty}^{+\infty} F(\mu) e^{j 2 \pi \mu t} d \mu f(t)=+F(μ)ej2πμtdμ
    其 中 , 2 π μ 代 表 频 率 , t 代 表 时 间 。 其中,2 \pi \mu代表频率,t代表时间。 2πμt傅里叶变换认为一个周期函数(信号)包含多个频率分量,任意函数(信号)f(t)可通过多个周期函数(基函数)相加而合成。从物理角度理解傅里叶变换是以一组特殊的函数(三角函数)为正交基,对原函数进行线性变换,物理意义便是原函数在各组基函数的投影

    1.2.二维傅里叶变换

    同样连续的傅里叶变换为:
    F ( μ , v ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ f ( t , z ) e − j 2 π ( μ t + v z ) d t d z F(\mu, v)=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(t, z) e^{-j 2 \pi(\mu t+v z)} d t d z F(μ,v)=++f(t,z)ej2π(μt+vz)dtdz
    它的反变换:
    f ( t , z ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ F ( μ , v ) e j 2 π ( μ t + v z ) d μ d v f(t, z)=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} F(\mu, v) e^{j 2 \pi(\mu t+v z)} d \mu d v f(t,z)=++F(μ,v)ej2π(μt+vz)dμdv
    傅里叶逆变换是将图像的频率分布函数变换为灰度分布函数傅里叶变换以前,图像(未压缩的位图)是由对在连续空间(现实空间)上的采样得到一系列点的集合,通常用一个二维矩阵表示空间上各点,记为z=f(x,y)。又因空间是三维的,图像是二维的,因此空间中物体在另一个维度上的关系就必须由梯度来表示,这样我们才能通过观察图像得知物体在三维空间中的对应关系。

    傅里叶频谱图上我们看到的明暗不一的亮点,其意义是指图像上某一点与邻域点差异的强弱,即梯度的大小,也即该点的频率的大小(可以这么理解,图像中的低频部分指低梯度的点,高频部分相反)。一般来讲,梯度大则该点的亮度强,否则该点亮度弱。这样通过观察傅里叶变换后的频谱图,也叫功率图,我们就可以直观地看出图像的能量分布:如果频谱图中暗的点数更多,那么实际图像是比较柔和的(因为各点与邻域差异都不大,梯度相对较小);反之,如果频谱图中亮的点数多,那么实际图像一定是尖锐的、边界分明且边界两边像素差异较大的。

    2.python×傅里叶级数

    理论过于抽象,还是实实在在的用代码折腾一下比较好理解。首先我们需要一个原函数:
    f ( x ) = { exp ⁡ ( x ) , − π ≤ x < 0 1 , 0 ≤ x < π f(x)=\left\{\begin{array}{lr} \exp (x), & -\pi \leq x<0 \\ 1, & 0 \leq x<\pi \end{array}\right. f(x)={exp(x),1,πx<00x<π
    函数图像如下:

    在这里插入图片描述

    可以用numpy和plt画出,代码如下:

    import os
    import cv2
    from pylab import *
    import matplotlib.pyplot as plt
    import numpy as np
    x = np.linspace(-3,pi)# 定义x轴
    cond = [True if (i>-pi and i<0) else False for i in x] #使用列表解析的方法
    y=1*(x>0)+np.exp(x)*cond # y的函数
    plt.plot(x,y)
    plt.show()
    

    接着算出该函数的a0,an,bn,最后自定义函数画出傅里叶级数变换,图像如下:

    在这里插入图片描述

    由上图看出,这图像就有正弦余弦内味了,周期函数都可以表示为不同频率的正弦和或者余弦和的形式。代码如下,建议做这种小实验使用jupyter,方便快捷

    import os
    import cv2
    from pylab import *
    import matplotlib.pyplot as plt
    import numpy as np
    x = mgrid[-10:10:0.2]
    n = arange(1, 1000)
    def fourier_transform():
        a0 = (1-exp(-pi))/pi+1
        s = a0 / 2
        for i in range(1, 100, 1):
            s0 = ((1-(-1)**i*exp(-pi))/(pi*(1+i**2))*cos(i*x)+1/pi*( (-i*(1-(-1)**i*exp(-pi)))/(1+i**2) + (1-(-1)**i)/i ) * sin(i*x))
            s = s + s0 
        plot(x,s,'orange',linewidth=0.6)
        title('fourier_transform')
        show() 
    fourier_transform()
    

    2.1.傅里叶变换后的频谱图

    使用傅里叶变换获得的频谱图。图像的主要成分是低频信息,它形成了图像的基本灰度等级,对图像结构的决定作用较小;高频信息形成了图像的边缘和细节,是在中频信息上对图像内容的进一步强化。
    在这里插入图片描述
    如上图所示,中间的亮点代表着图像的低频,其余周围则表示图像的高频信息。若不将低频信息移至中间,则低频信息过于分散,频谱未中心化如下图所示:
    在这里插入图片描述
    代码如下所示:

    import cv2
    import numpy as np
    import matplotlib.pyplot as plt
    # 读取图片
    img = cv2.imread(' ', 0)
    
    # 进行float32形式转换
    float32_img = np.float32(img)
    
    # 使用cv2.dft进行傅里叶变化
    dft_img = cv2.dft(float32_img, flags=cv2.DFT_COMPLEX_OUTPUT)
    
    # 使用np.fft.shiftfft()将变化后的图像的低频转移到中心位置
    dft_img_ce = np.fft.fftshift(dft_img)
    
    # 使用cv2.magnitude将实部和虚部转换为实部,乘以20是为了使得结果更大
    # img_dft = 20 * np.log(cv2.magnitude(dft_img_ce[:, :, 0], dft_img_ce[:, :, 1]))
    # 频谱未中心化
    img_dft = 20 * np.log(cv2.magnitude(dft_img[:, :, 0], dft_img[:, :, 1]))
    # 进行画图操作
    plt.subplot(121) # 我是第一行第二列的第一个
    plt.imshow(img, cmap='gray')
    plt.subplot(122)
    plt.imshow(img_dft, cmap='gray')
    plt.show()
    

    3.频率域滤波

    3.1.低频与高频

    低频对应图像内变换缓慢的灰度分两。例如在一幅大草原的图像中,低频对应着广袤的颜色趋于一致的草原。总结就是低频是图像中平坦的,灰度变化不大的点,图像中的大部分区域

    高频对应图像内变化越来越快的灰度分量,是由灰度的尖锐过度造成的,例如,还是大草原那幅图像,若图像中有一只狮子,则狮子的边缘信息就是高频。总结就是高频是图像中灰度变化剧烈的点,一般是图像轮廓或者是噪声

    根据图像的高频与低频的特征,我们可以设计相应的高通与低通滤波器。通过低频的滤波器叫做低通滤波器,通过高频的滤波器叫做高通滤波器。高通滤波可以检测图像中尖锐、变化明显的地方;低通滤波可以让图像变得光滑,滤除图像中的噪声

    3.2.频率域滤波步骤

    为了方便下述低通高通滤波器的实验,必须先明确频率域滤波的基本步骤。

    1. 首先将图片读取成灰度图且转化为np.float32()
    2. 使用cv2.dft函数进行傅里叶变化
    3. 使用np.fft.fftshift函数将低频转移到图像中心
    4. 定义掩模:生成的掩模中间为1周围为0
    5. 将掩模与傅里叶变化后图像相乘,保留中间部分
    6. 使用np.fft.ifftshift将低频移动到原来的位置
    7. 使用cv2.idft进行傅里叶的反变化
    8. 使用cv2.magnitude转化为空间域内
    9. 使用plt进行绘图展示

    由上述叙述可知,滤波的关键取决于滤波函数,也称为滤波器,因为它在滤波中抑制或除去了频谱中的某些分量,而保留其他的一些频率不受影响,从而达到滤波的目的。而下述的 低通滤波器和高通滤波器的不同点,便在于定义不同的滤波,即步骤中的掩模,不同的滤波将产生不同的效果。

    3.3.低通滤波器

    理想的低通滤波器模板为:
    H ( u , v ) = { 1 , D ( u , v ) ≤ D 0 0 , D ( u , v ) > D 0 H(u, v)=\left\{\begin{array}{ll} 1, & D(u, v) \leq D_{0} \\ 0, & D(u, v)>D_{0} \end{array}\right. H(u,v)={1,0,D(u,v)D0D(u,v)>D0
    其中,D0表示通带半径,D(u,v)是到频谱中心的距离(欧式距离),计算公式如下:
    D ( u , v ) = ( u − M / 2 ) 2 + ( v − N / 2 ) 2 \mathrm{D}(u, v)=\sqrt{(u-M / 2)^{2}+(v-N / 2)^{2}} D(u,v)=(uM/2)2+(vN/2)2
    M和N表示频谱图像的大小,(M/2,N/2)即为频谱中心。低通滤波器代码如下所示:

    # 低通滤波器
    import cv2
    import numpy as np
    from matplotlib import pyplot as plt
    # 第一步读入图片
    img = cv2.imread(' ', 0)# filepath
    # 使用cv2.dft进行傅里叶变化
    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    # 使用np.fft.fftshift将低频转移到图像中心
    dft_center = np.fft.fftshift(dft)
    # 定义掩模:生成的掩模中间为1周围为0
    crow, ccol = int(img.shape[0] / 2), int(img.shape[1] / 2) # 求得图像的中心点位置
    mask = np.zeros((img.shape[0], img.shape[1], 2), np.uint8)
    mask[crow-30:crow+30, ccol-30:ccol+30] = 1
    
    # 将掩模与傅里叶变化后图像相乘,保留中间部分
    mask_img = dft_center * mask
    
    # 使用np.fft.ifftshift(将低频移动到原来的位置
    img_idf = np.fft.ifftshift(mask_img)
    
    # 使用cv2.idft进行傅里叶的反变化
    img_idf = cv2.idft(img_idf)
    
    # 使用cv2.magnitude转化为空间域内
    img_idf = cv2.magnitude(img_idf[:, :, 0], img_idf[:, :, 1])
    
    # 进行绘图操作
    plt.subplot(121)
    plt.imshow(img, cmap='gray')
    plt.subplot(122)
    plt.imshow(img_idf, cmap='gray')
    plt.show()
    

    在这里插入图片描述

    由实验可以看出,图片经过低通滤波器,使低频通过而使高频衰减 ,降低了图像变换的幅度。平滑且模糊掉了尖锐和边缘部分。保留了图像的整体概貌(低频部分),去掉了图像的细节(高频部分)

    3.4.高通滤波器

    和低通滤波器相反,理想的高通滤波器为:
    H ( u , v ) = { 0 , D ( u , v ) ≤ D 0 ) 1 , D ( u , v ) > D 0 ) H(u, v)=\left\{\begin{array}{ll} 0, & \left.\mathrm{D}(\mathrm{u}, \mathrm{v}) \leq D_{0}\right) \\ 1, & \left.\mathrm{D}(\mathrm{u}, \mathrm{v})>D_{0}\right) \end{array}\right. H(u,v)={0,1,D(u,v)D0)D(u,v)>D0)
    其中,D0表示通带半径,D(u,v)是到频谱中心的距离(欧式距离),计算公式如下:
    D ( u , v ) = ( u − M / 2 ) 2 + ( v − N / 2 ) 2 \mathrm{D}(u, v)=\sqrt{(u-M / 2)^{2}+(v-N / 2)^{2}} D(u,v)=(uM/2)2+(vN/2)2

    常见的巴特沃斯高通滤波器(BHPF)公式如下:
    H ( u , v ) = 1 1 + [ D 0 / D ( u , v ) ] 2 n H(u, v)=\frac{1}{1+\left[D_{0} / D(u, v)\right]^{2 n}} H(u,v)=1+[D0/D(u,v)]2n1
    还有高斯高通滤波器(GHPF):
    H ( u , v ) = 1 − e − D 2 ( u , v ) / 2 σ 2 , σ = D 0 H(u, v)=1-e^{-D^{2}(u, v) / 2 \sigma^{2}}, \sigma=D_{0} H(u,v)=1eD2(u,v)/2σ2,σ=D0

    # 实现高通滤波器
    import cv2
    import numpy as np
    from matplotlib import pyplot as plt
    
    #读取图像
    img = cv2.imread(' ', 0) # filepath
    #傅里叶变换
    dft = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)
    fshift = np.fft.fftshift(dft)
    #设置高通滤波器
    rows, cols = img.shape
    crow,ccol = int(rows/2), int(cols/2) #中心位置
    mask = np.ones((rows, cols, 2), np.uint8) # 低通滤波器这是np.zeros
    mask[crow-30:crow+30, ccol-30:ccol+30] = 0 # 低通滤波器这就是1
    #掩膜图像和频谱图像乘积
    f = fshift * mask
    #傅里叶逆变换
    ishift = np.fft.ifftshift(f)
    iimg = cv2.idft(ishift)
    res = cv2.magnitude(iimg[:,:,0], iimg[:,:,1])
    #显示原始图像和高通滤波处理图像
    plt.subplot(121), plt.imshow(img, 'gray')
    plt.subplot(122), plt.imshow(res, 'gray')
    plt.show()
    

    实现效果如下图所示:

    在这里插入图片描述

    由上述实验可以看出:高通滤波器使高频通过而使低频衰减,被高通滤波的图像比原始图像少灰度级的平滑过渡而突出边缘等细节部分,保留了图像的细节,模糊了概貌。也就是边缘信息更加突出,例如轮廓

    3.5.低通与高通滤波器实验总结

    由上述实验可看出,低通滤波器的滤波定义为:

    crow, ccol = int(img.shape[0] / 2), int(img.shape[1] / 2) # 求得图像的中心点位置
    mask = np.zeros((img.shape[0], img.shape[1], 2), np.uint8)
    mask[crow-30:crow+30, ccol-30:ccol+30] = 1
    

    高通滤波器的滤波为:

    crow, ccol = int(img.shape[0] / 2), int(img.shape[1] / 2) # 求得图像的中心点位置
    mask = np.ones((img.shape[0], img.shape[1], 2), np.uint8) # 低通滤波器这是np.zeros
    mask[crow-30:crow+30, ccol-30:ccol+30] = 0 # 低通滤波器这就是1
    

    相同点在于定义的掩模大小和图片大小相同,不同点在于低通滤波器中间为1,那么经过傅里叶变换的图像和掩模相乘之后,便只保留中间部分;相反,高通滤波器中间为0,经过傅里叶变换的图像和掩模相乘之后,便只保留边缘部分

    展开全文
  • 基于MATLAB图像处理的频率域滤波分析及应用
  • 本文档包含数字图像处理(冈萨雷斯第三版)中的图像频率域滤波(傅里叶变换、巴特沃斯低通滤波高通滤波等)实验,包含整个实验过程和原理解释,以及详细的执行代码。代码复制后可在matlab中直接运行。
  • 一)空间域滤波与频率域滤波  1)空间域滤波  空间域滤波是指在图像空间中借助模板对图像领域进行操作,处理图像每一个像素值。主要分为线性滤波和非线性滤波两类,根据功能可分为平滑滤波器和锐化滤波器。...
  • 基于MATLAB图像处理的频率域滤波分析及应用.pdf
  • 频率域滤波的MATLAB设计与实现课程设计报告书.doc
  • 频率域滤波去除周期性噪声

    千次阅读 2020-12-26 03:25:40
    频率域滤波去除周期性噪声(波纹状噪声)matlab实现
  • 频率域滤波(1)

    2019-10-02 11:44:15
    频率域滤波实际上是将图像进行傅里叶变换,然后在变换域进行处理,然后进行傅里叶反变换转换回空间域,原理是用傅里叶变换表示的函数特征完全可以通过傅里叶反变换来重建,而且不会丢失任何信息(因为任何周期或非...
  • 7.频率域滤波基础

    2019-08-23 07:03:09
    目录 一 原理 ...三 频率域滤波的计算步骤 步骤1-2解释 四 示例 五 频率域滤波与空间域滤波的关系 1.给出一个频率域的滤波器,如果得到其空间域的滤波器? 2.到底用哪个? 怎么用? 3.举...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 8,196
精华内容 3,278
关键字:

频率域滤波