数字图像学笔记——10. 频域与傅里叶分析方法

水深无声 2022-11-07 13:30 282阅读 0赞

> 频域滤波技术，目前主要使用的有两种类型，一种是傅立叶变换技术，还有一种是小波分析。基本逻辑就是把原始信号映射到频率空间中，使得在时域空间无法处理的信号，得以在另外一种空间体系下能够被较有效的处理。

> 目前，我已经完成了全部的傅立叶分析的科普性的文章，所以在这个系列里我就不再赘述这方面的基础知识，我会在这下面贴出这些文章的链接，你可以根据自己的需要去查阅对应的章节，如果想直接了解如何实现一个一维傅立叶FFT与iFFT算法的朋友，直接从6、7、8三个章节看就行了。

[《漫谈傅里叶1——从无穷级数到傅里叶》][1]  
[《漫谈傅里叶2——公式推导、三角函数正交性》][2]  
[《漫谈傅里叶3——收敛性、非周期函数的推广应用》][3]  
[《漫谈傅里叶4——全时傅里叶的缺点与短时傅里叶》][4]  
[《漫谈傅里叶5——卷积与短时傅里叶的缺点》][5]  
[《漫谈傅里叶6——采样与1D初步实现》][6_1D]  
[《漫谈傅里叶7——带有相位与幅值的1D实现》][7_1D]  
[《漫谈傅里叶8——傅里叶逆计算的实现》][8]

而关于小波分析这部分，我打算近期更新一些关于量化金融方面的文章后再回头来写，你可以关注我的博客，或者新浪微博——打码的阿通。

### 文章目录 ###

*  用OpenCV实现的图片DFT示例
 *   *  正向 DFT 或 FFT 变换后的结果
     *  逆向DFT 或从频率转回图像的结果
 *  理想滤波器 (Ideal Filter)
 *  Butterworth 滤波器

# 用OpenCV实现的图片DFT示例 #

尽管我知道Python有很多方便的FFT库，但是为了方便学习C/C++的朋友，我尽量用原生的OpenCV库做完了这个用例。理论的那一套在前面的一维FFT里已经讲过，所以这里不赘述，直接上完整的代码好了。

import cv2
    
    from Common import *
    
    
    def load_img_grayscale(file: str):
        image = cv2.imread(file, cv2.IMREAD_GRAYSCALE)
        return image
    
    
    def shift_spectrum(gamma_matrix):
    	# 和C++语言的OpenCV库中获取矩阵大小方式类似
    	# 可以用如下方法获得
    	# int rows = gamma_matrix.rows
    	# int cols = gamma_matrix.cols
        rows, cols = gamma_matrix.shape
        cr = int(rows / 2)
        cc = int(cols / 2)
    
    	# 计算四个象限子矩阵
    	# 可以使用C/C++中通常使用 Rect 来获取原矩阵的四个象限范围
    	# (0,0)----------------------------->cols
    	# Q2 | Q1 |
    	# (0r,0c) | | 
    	# |(cr,cc) |
    	# ---------------|------------------- |
    	# Q3 | Q4 |
    	# | |
    	# | |
    	# (rows,cols) 
    	
        quadrant_1 = gamma_matrix[0:cr, cc:cols]
        quadrant_2 = gamma_matrix[0:cr, 0:cc]
        quadrant_3 = gamma_matrix[cr:rows, 0:cc]
        quadrant_4 = gamma_matrix[cr:rows, cc:cols]
    
        temp = quadrant_1.copy()
        quadrant_1 = quadrant_3
        quadrant_3 = temp
    
        temp = quadrant_2.copy()
        quadrant_2 = quadrant_4
        quadrant_4 = temp
    
        # 将修改好的数值写入到output矩阵中
        output = np.zeros(gamma_matrix.shape)
        output[0:cr, cc:cols] = quadrant_1
        output[0:cr, 0:cc] = quadrant_2
        output[cr:rows, 0:cc] = quadrant_3
        output[cr:rows, cc:cols] = quadrant_4
    
        return output
    
    
    def gen_spectrum(real_matrix, imag_matrix):
        # 创建一个gamma变化后的矩阵
        gamma_matrix = np.zeros(real_matrix.shape, dtype=np.float32)
        cv2.magnitude(real_matrix, imag_matrix, gamma_matrix)
        cv2.log(gamma_matrix + 1, gamma_matrix)
    
        # 归一化：normalize
        output = np.zeros(real_matrix.shape, dtype=np.float32)
        cv2.normalize(gamma_matrix, output, 0, 1, cv2.NORM_MINMAX)
    
        return output
    
    
    def do_fft():
        # 获取原始图像的灰度图
        img = load_img_grayscale("Data/Illustrations/1.jpg")
    
        # 获取图像的长宽
        # C++代码里，需要用 img.rows, img.cols来获取对应的数据
        rows, cols = img.shape
    
        # 对长宽进行优化，DFT算法需要输入的数组长度为2^n
        m = cv2.getOptimalDFTSize(rows)
        n = cv2.getOptimalDFTSize(cols)
    
        # 按照要求，对数据进行填充，不足的部分填充0
        # C++中使用
        # cv2.copyMakeBorder(src, dst, 0, m - rows, 0, n - cols, cv2.BORDER_CONSTANT)
        fft2 = np.zeros((m, n, 2), np.float32)
        fft2[:rows, :cols, 0] = img
    
        # 为DFT创建实数部和虚数部
        # C++对应的方法：
        # Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
        # Mat complex;
        # cv2::merge(planes, 2, complex):
        # dft(complexImg, complexImg, cv2.DFT_COMPLEX_OUTPUT)
    
        # 使用DFT算法，进行快速傅立叶变换
        cv2.dft(fft2, fft2, cv2.DFT_COMPLEX_OUTPUT)
    
        # 分离复数矩阵
        # C/C++中，用的是 split(fft2, real_mat, imag_mat) 进行分离
        real_mat, imag_mat = cv2.split(fft2)
    
        # 生成频率数据，由于一部分数据属于不可见数据，所以需要进行gamma变换
        # 并执行归一化操作
        gamma = gen_spectrum(real_mat, imag_mat)
    
        # 为gamma图执行中央转换
        # 方法是1，3象限对调，2，4象限对调
        shift = shift_spectrum(gamma_spectrum)
    
        # 把原始图像、FFT频率图、转化后的频率图全部显示出来
        plt = PltImageCache()
        plt.add(img, "origin")
        plt.add(gamma, "dft computed")
        plt.add(shift, "dft shifted")
        plt.plots(1, 3)
    
        # 返回复数矩阵供下一步操作
        return fft2
    
    
    def do_ifft(fft2):
    
        # 获取原始图像的灰度图
        img = load_img_grayscale("Data/Illustrations/1.jpg")
    
        ifft2 = np.zeros(fft2.shape[:2], np.float32)
        cv2.dft(fft2, ifft2, cv2.DFT_REAL_OUTPUT | cv2.DFT_INVERSE | cv2.DFT_SCALE)
    
        # 把原始图像、FFT频率图、转化后的频率图全部显示出来
        plt = PltImageCache()
        plt.add(img, "original")
        plt.add(ifft2, "idft")
        plt.plots(1, 2)
    
    
    if __name__ == "__main__":
        dft_matrix = do_fft()
        do_ifft(dft_matrix)

输出的结果就是这样的了：

## 正向 DFT 或 FFT 变换后的结果 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center]

## 逆向DFT 或从频率转回图像的结果 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center 1]

这里你会发现转换回来的照片多了一圈黑色的边框，那个就是为了进行DFT运算，而调整的新图片大小，长宽大小都是  2 n 2^n 2n。

通常，在做完FFT运算后，它输出的结果是对称的，对于2维的也是这样，不做任何处理的频率图，频率范围从四边角到中央逐渐递增，所以为了方便观察，会把四个象限的子矩阵做一下调整。

所以最终的图片也就成了从中央一个点向四周呈现放射状的形式。因此频率范围的调整就变成了计算圆半径R的问题。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center 2]  
通常，我们为了过滤图片中的一些信息，就会使用一些滤波器，下面介绍一些比较常用的频域滤波器。

# 理想滤波器 (Ideal Filter) #

类似于二极管的作用，只放走设定域值内的频率，而域值外的频率一律设置为0

f ( x ) = \{ g ( x ) a < x ≤ b 0 e l s e f(x) = \\left\\\{\\begin\{matrix\} g(x) & a < x \\leq b \\\\ 0 & else \\\\ \\end\{matrix\}\\right. f(x)=\{ g(x)0a<x≤belse

而在设代码结构上与椒盐噪音的方式是相似的，其实现代码大概是这样的：

def generate_ideal_mask(rows, cols, alpha, beta=0):
        """ :param rows: rows of spectrum :param cols: columns of spectrum :param alpha: the frequency upper limit :param beta: the frequency lower limit :return: """
    
        import numpy as np
    
        output = np.zeros((rows, cols), np.float32)
    
        # center coordinate
        cr = rows / 2
        cc = cols / 2
    
        for r in range(rows):
            for c in range(cols):
                distance = np.sqrt((r - cr) ** 2 + (c - cc) ** 2)
                if alpha > distance >= beta:
                    output[r, c] = 1.
                else:
                    output[r, c] = 0.
    
        return output

而输出结果如下：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center 3]

# Butterworth 滤波器 #

其数学解析式如下：

f ( x ) = 1 1 + \[ D ( u , v ) D o \] 2 n f(x) = \\frac\{1\}\{1 +\[\\frac\{D(u,v)\}\{D\_o\}\]^\{2n\}\} f(x)=1\+\[DoD(u,v)\]2n1

代码的实现方式大体如下：

def generate_butterworth_mask(rows, cols, n, d, flip=False):
        """ :param rows: rows of spectrum :param cols: columns of spectrum :param n: the filter adjustment factor :param d: the d0 :return: """
    
        import numpy as np
    
        output = np.zeros((rows, cols), np.float32)
    
        # center coordinate
        cr = rows / 2
        cc = cols / 2
    
        for r in range(rows):
            for c in range(cols):
                distance = np.sqrt((r - cr) ** 2 + (c - cc) ** 2)
    
                frac = (distance / d) ** (2 * n)
                output[r, c] = 1 / (1 + frac)
    
        if flip:
            output = 1 - output
    
        return output

得到的效果大概是这样的：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center 4]

[1]: https://seagochen.blog.csdn.net/article/details/80840211
[2]: https://seagochen.blog.csdn.net/article/details/80896471
[3]: https://seagochen.blog.csdn.net/article/details/80914127
[4]: https://seagochen.blog.csdn.net/article/details/81004560
[5]: https://seagochen.blog.csdn.net/article/details/81029266
[6_1D]: https://seagochen.blog.csdn.net/article/details/113789264
[7_1D]: https://seagochen.blog.csdn.net/article/details/114127894
[8]: https://seagochen.blog.csdn.net/article/details/114275642
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center]: /images/20221023/45dcd3875dcb41a98b6292835da5675f.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221023/3b405510ec644e9b8b5c3b6f990aee4d.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center 2]: /images/20221023/7836e16568b248dcb243ce96f246b840.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center 3]: /images/20221023/74f74b3454db47879e82da4da1c96275.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3BvaXNvbmNocnk_size_16_color_FFFFFF_t_70_pic_center 4]: /images/20221023/01a138ad6464435793071aba4caf82f7.png