数字图像学笔记——12. 图像退化与复原（频域噪音函数对图像的影响、以及什么是陷滤波器）

ゝ一世哀愁。 2022-08-28 12:59 139阅读 0赞

> 聊完了空间噪音对函数的一般影响后，我们需要来聊一聊由于频域上的噪音对图像的影响了。

### 文章目录 ###

*  频率噪音产生的原因
 *  频域噪音在图像上表现的特点
 *   *  通过DFT函数获取图像的频谱图
     *  有频域噪音信号的图像的特点
 *  频域滤波器——陷滤波器（Notch Filter）
 *  实现代码
 *  相关示例代码

# 频率噪音产生的原因 #

当我们的视觉采集设备从传统的化学方法（卤化银显影法）变成了后来的数码（CMOS）设备后，由于CMOS作为一种场效应晶体管的原因，所以会受到温度、湿度、磁场环境的影响。所以在我们的采集设备都变成了数码设备后，由于电子产品的通病，所以才开始有了频域上的噪音。

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

# 频域噪音在图像上表现的特点 #

## 通过DFT函数获取图像的频谱图 ##

如果通过numpy实现这个功能会非常简单

def frequency_noise_analysis(filepath: str):
        img = load_image_gray(filepath)
    
        # convert byte to float
        dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    
        # use NumPy to rapidly shift DFT diagram and prepare for display FFT diagram
        dft_shift = np.fft.fftshift(dft)
        result = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
    
        return [img, result]

当然你也可以通过cv自带的函数来获取频谱图

def shift_spectrum(gamma_matrix):
        rows, cols = gamma_matrix.shape
        cr = int(rows / 2)
        cc = int(cols / 2)
    
        quadrant_1 = gamma_matrix[cr:rows, 0:cc]
        quadrant_2 = gamma_matrix[0:cr, 0:cc]
        quadrant_3 = gamma_matrix[0:cr, cc:cols]
        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
    
        # yield the updates value back
        output = np.zeros(gamma_matrix.shape)
        output[cr:rows, 0:cc] = quadrant_1
        output[0:cr, 0:cc] = quadrant_2
        output[0:cr, cc:cols] = quadrant_3
        output[cr:rows, cc:cols] = quadrant_4
    
        return output
    
    def do_fft(imgfile: str):
        # 获取原始图像的灰度图
        img = cv2.imread(imgfile, cv2.IMREAD_GRAYSCALE)
    
        # 获取图像的长宽
        # 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)
    
        # 分离复数矩阵
        real_mat, imag_mat = cv2.split(fft2)
    
        # 生成频率数据，由于一部分数据属于不可见数据，所以需要进行gamma变换
        # 并执行归一化操作
        gamma = gen_spectrum(real_mat, imag_mat)
    
        # 为gamma图执行中央转换
        # 方法是1，3象限对调，2，4象限对调
        shift = shift_spectrum(gamma)
    
        # 把原始图像、FFT频率图、转化后的频率图全部显示出来
        plt = DiagramPlotter()
        plt.append_image(img, "origin")
        plt.append_image(gamma, "dft computed")
        plt.append_image(shift, "dft shifted")
        plt.show(1, 3)
    
        # 返回复数矩阵供下一步操作
        return fft2

如果你主打C、Java语言做图像处理方面的工作，相信以上代码理解以后可以很容易改成C、Java语言的。

## 有频域噪音信号的图像的特点 ##

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

一个在频域上有噪音的图像，会表现出明显的周期噪音信号的特点，如果把图片做FFT分析后，会发现围绕着中心点周边出现数个集中的疑似噪点信号，请记住对于二维DFT来说，通常中间附近都是低频信号，而越往边缘方向，越是高频信号，低频信号一般包含照片细节信息，而照片边缘信息则通常处于高频区。

# 频域滤波器——陷滤波器（Notch Filter） #

所以，为了消除频域信号上的噪音，我们通常会想到直接把频域信号上的噪点用某种方法滤除掉，所以自然就会想到利用类似掩码的方式滤除频域噪音，或者说用类似挖西瓜瓤的方式去除那些个噪点，而这在信号学领域有个专门的名称，就是——陷滤波器。

![在这里插入图片描述][watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBA5omT56CB55qE6Zi_6YCa_size_20_color_FFFFFF_t_70_g_se_x_16_pic_center 1]  
它就像个勺子，专门用来挖掉频谱图上你怀疑的异常高亮点。

由于使用滤波器，大多数情况下都是手动工作（顺便说一下，我最近正在研究怎么用AI自适应噪音信号滤除，有进展后我会写一篇博文介绍一下相关工作）。根据我们多年的工作经验，对于上图的噪音，极有可能是除中心点外那8个点异常亮的地方。

所以一个实现策略，就是直接做一个类似MASK的东西，把黑色的地方全部删除掉，然后再还原回图像。

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

接下来我们就来看看怎么实现吧。

# 实现代码 #

这里面我们最主要的是要实现这个代码

def use_notch_filter(dft, x, y, radius):
        thetas = np.linspace(-np.pi, np.pi, 100)
        
        for r in range(-radius, radius, 1):
            for theta in thetas:
                xpt = r * np.cos(theta) + x
                ypt = r * np.sin(theta) + y
    
                # round number to integer
                xpt = int(round(xpt))
                ypt = int(round(ypt))
    
                # assign value
                dft[xpt, ypt] = 0
    
        return dft

它的作用是在我们指定的位置上生成一个圆形的区域，然后所有在这个区域内的数据全部被置零。然后在我们的图片进过FFT后，手动设置一下坐标和大小，就可以了。

def frequency_noise_analysis(filepath: str):
        # load image from file
        img = load_image_gray(filepath)
    
        # convert byte to float
        dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
    
        # use NumPy to rapidly shift DFT diagram and prepare for display FFT diagram
        dft_shift = np.fft.fftshift(dft)
    
        dft_filter = use_notch_filter(dft_shift.copy(), 40, 55, 5)      # filter noise pt 1
        dft_filter = use_notch_filter(dft_filter, 80, 55, 5)      # filter noise pt 2
        dft_filter = use_notch_filter(dft_filter, 165, 55, 5)     # filter noise pt 3
        dft_filter = use_notch_filter(dft_filter, 205, 55, 5)     # filter noise pt 4
    
        dft_filter = use_notch_filter(dft_filter, 40, 110, 5)      # filter noise pt 1
        dft_filter = use_notch_filter(dft_filter, 80, 110, 5)      # filter noise pt 2
        dft_filter = use_notch_filter(dft_filter, 165, 110, 5)     # filter noise pt 3
        dft_filter = use_notch_filter(dft_filter, 205, 110, 5)     # filter noise pt 4
    
        # convert dft image back
        img_back = np.fft.fftshift(dft_filter)
        img_back = cv2.idft(img_back, flags=cv2.DFT_COMPLEX_INPUT | cv2.DFT_SCALE | cv2.DFT_REAL_OUTPUT)
    
        # prepare image
        dft_img = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))
        filter_img = 20 * np.log(cv2.magnitude(dft_filter[:, :, 0], dft_filter[:, :, 1]))
    
        return [img, dft_img, filter_img, img_back]

最终的实现效果：

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

是不是好多了呢。

# 相关示例代码 #

经常有朋友在后台私信我要示例代码，这次我就把代码贴在这里，相关的内容已经更新并推送到了[Github(Python 源码)][Github_Python]，有需要的朋友就自己下载吧。

如果你觉得这个内容对你有帮助，点个赞再走呗~~

[watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBA5omT56CB55qE6Zi_6YCa_size_19_color_FFFFFF_t_70_g_se_x_16_pic_center]: /images/20220828/03a034c182a8486fa03518b40b44d00b.png
[watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBA5omT56CB55qE6Zi_6YCa_size_20_color_FFFFFF_t_70_g_se_x_16_pic_center]: /images/20220828/7fb11be7808a460893e1a711d084a76d.png
[watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBA5omT56CB55qE6Zi_6YCa_size_20_color_FFFFFF_t_70_g_se_x_16_pic_center 1]: /images/20220828/eb24de07496b4e62a9d040a30c261cc7.png
[watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBA5omT56CB55qE6Zi_6YCa_size_19_color_FFFFFF_t_70_g_se_x_16_pic_center 1]: /images/20220828/e35f300d7bba407eb306d325d6f0e22e.png
[watermark_type_ZHJvaWRzYW5zZmFsbGJhY2s_shadow_50_text_Q1NETiBA5omT56CB55qE6Zi_6YCa_size_20_color_FFFFFF_t_70_g_se_x_16_pic_center 2]: /images/20220828/e1bf67a1b74348cea58461bc41e73359.png
[Github_Python]: https://github.com/seagochen/AlgorithmLearning/blob/master/Cv/Degradation/FrequencyNoise.py