数字图像处理知识点梳理——第四章 频率域滤波

理论部分

本章和第三章一样，主要应用于图像增强，内容大致分为以下这几个部分：
一、傅里叶变换及其相关知识点
在维度上分为一维傅里叶变换和二维傅里叶变换
在连续性上分为连续傅里叶变换和离散傅里叶变换
二、频率域滤波器
根据频率的高低分为低通滤波器和高通滤波器

一、傅里叶变换及其相关知识点

傅里叶技术为研究和实现主要的图像处理算法提供了有意义且实用的方法。在本章的学习中，我们首先要了解到几个相关的名词：

冲激：或者称为分布或者广义函数，冲激可以是连续的或离散的，连续的用$\int$，离散的用$\sum$，以下只列出了连续情况：
一维情况下：对于任意点 $t_0$ 的冲激，表示为 $\delta (t-t{}_0)$
二维情况下：对于位于坐标$(t_0,z_0)$处的冲激，表示为 $\delta (t-t_0,z-z_0)$

取样特性：连续的用$\int$，离散的用$\sum$，以下只列出了连续情况：
一维情况下：对于任意点 $t_0$ 的取样特性表示为$${\int }_{-\infty }^{+\infty }f(t)\delta (t-t{}_0)dt=f(t_0)$$
二维情况下：对于位于坐标$(t_0,z_0)$处的取样特性为$${\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }f(t,z)\delta (t-t_0,z-z_0)dtdz=f(t_0,z_0)$$
冲激串：就是对冲激进行了累加
一维情况下：表示为$$s_{\Delta T}(t)=\sum _{n=-\infty }^{\infty }\delta (t-n\Delta T)$$，周期为$\Delta T$的冲激串的傅里叶变换还是冲激串，其周期为$1/\Delta T$，取样间隔要满足$\Delta T<1/2\mu$
二维情况下：$$s_{\Delta T\Delta Z}(t,z)=\sum _{m=-\infty }^{\infty }\sum _{n=-\infty }^{\infty }\delta (t-m\Delta T,z-n\Delta Z)$$,取样间隔要满足$\Delta T<1/2{\mu }_{max}且\Delta Z<1/2v_{max}$
欧拉公式：$$e^{j\theta }=cos\theta +jsin\theta$$
sinc函数：$$sinc(m)=\frac{sin(\pi m)}{(\pi m)}$$
卷积：在第三章的时候接触过
连续情况下记为：$f(t)★h(t)={\int }_{-\infty }^{+\infty }f(\tau )h(t-\tau )d\tau$；
离散情况下记为：$f(x)★h(x)={\sum }_{m=0}^{M-1}f(m)h(x-m)$

对$f(t)★h(t)$进行傅里叶变换，可得$\zeta [f(t)★h(t)]=H(\mu )F(\mu )$；
$f(t)★h(t)$和$H(\mu )F(\mu )$是傅里叶变换对，可以写为：$$f(t)★h(t)⟺H(\mu )F(\mu )$$（一维情况下）$$f(x,y)★h(x,y)⟺F(u,v)H(u,v)$$（二维情况下）

卷积定理的另一半：$$f(t)h(t)⟺H(\mu )★F(\mu )$$（一维情况下）
$$f(x,y)h(x,y)⟺F(u,v)★H(u,v)$$（二维情况下）

由此可以说明频率域的卷积类似于空间域的乘积，两者分别与傅里叶正、反变换相联系。
学第三章的时候，不是很清楚相关和卷积是用来做什么的，现在了解了，卷积定理是频率域滤波的基础。

取样定理：如果以超过函数最高频率的两倍的取样率来获得样本，连续的带限函数可以完全的从他的样本集来恢复，这就是取样定理。取样定理是数字信号处理理论的基石。
混淆：如果一个带限函数用低于其最高频率的两倍取样率取样将会发生周期重叠，不管使用什么滤波器，都不可能分离出变换的一个单周期，这就是频率混淆，简称为混淆。
在图像中，存在两种主要的混淆现象：空间混淆和时间混淆。空间混淆由欠取样造成，时间混淆与图像序列中图像间的时间间隔有关。
抗混淆：在取样过的函数中，混淆总是存在的，在实践中，我们可以通过平滑输入函数减少高频分量的方法来降低混淆的影响，这就是所谓的抗混淆，由于抗混淆是一个取样问题，所以它必须在被取样之前完成。
傅里叶变换的频谱、功率谱、相角：
一维情况下：
频谱：$|F(u)|=[R^2(u)+I^2(u){]}^{\frac{1}{2}}$
功率谱：$E(u)=|F(u){|}^2=R^2(u)+I^2(u)$
相角：$\varphi (u)=arctan\frac{I(u)}{R(u)}$
二维情况下：
频谱：$|F(u,v)|=[R^2(u,v)+I^2(u,v){]}^{\frac{1}{2}}$
功率谱：$E(u,v)=R^2(u,v)+I^2(u,v)$
相角：$\varphi (u,v)=arctan\frac{I(u,v)}{R(u,v)}$

（一）一维情况下

连续的傅里叶变换

傅里叶变换：$F(\mu )={\int }_{-\infty }^{+\infty }f(t)e^{-j2\pi \mu t}dt$
傅里叶反变换：$f(t)={\int }_{-\infty }^{+\infty }F(\mu )e^{j2\pi \mu t}d\mu$
傅里叶变换和傅里叶反变换合起来称为傅里叶变换对，以下离散情况和二维连续、离散情况也是如此称呼，不再赘述。

用欧拉公式带入$F(\mu )$，可以把式子变换为$$F(\mu )={\int }_{-\infty }^{+\infty }f(t)[cos(2\pi \mu t)-jsin(2\pi \mu t)]dt$$，类似的可以自己得到$f(t)$，其中，频率变量$\mu$的单位取决于$t$的单位

$但是在计算机处理之前，连续函数需要转换为离散值序列$，所以对连续的函数进行取样；
取样后的函数为$$\tilde{f}(t)=f(t)s_{\Delta T}(t)=\sum _{n=-\infty }^{\infty }f(t)\delta (t-n\Delta T)$$

取样函数的傅里叶变换为：$$\tilde{F}(t)=F(\mu )★S(\mu )=\frac{1}{\Delta T}\sum _{n=-\infty }^{\infty }F(\mu -\frac{n}{\Delta T})$$
$\tilde{F}(t)$的拷贝间隔由$1/\Delta T$决定，值过大时是过取样，会提供更清晰的间隔，值正好时是临界取样，值过小时是欠取样，导致周期融合。

离散的傅里叶变换（DFT）

对一个周期取样是DFT的基础，与连续的傅里叶变换做对比，离散的情况如下：
傅里叶变换：$F_m={\sum }_{n=0}^{M-1}{f}_n{e}^{-j2\pi mn/M}$
傅里叶反变换（IDFT)：$f_n=\frac{1}{M}{\sum }_{m=0}^{M-1}{F}_m{e}^{j2\pi mn/M}$
$m=0,1,2,...,M-1$

离散频率域中的相应间隔$\Delta u=\frac{1}{M\Delta T}=\frac{1}{T}$
DFT的频率分辨率$\Delta u$取决于连续函数$f(t)$被取样的持续时间$T$，并且DFT跨越的频率范围取决于取样间隔$\Delta T$

（二）二维情况下

连续的傅里叶变换

傅里叶变换：$F(\mu ,v)={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }f(t,z)e^{-j2\pi (\mu t+vz)}dtdz$
傅里叶反变换：$f(t,z)={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }F(\mu ,v)e^{j2\pi (\mu t+vz)}d\mu dv$
$\mu ,v$是频率变量，$\mu ,v$的域定义了连续频率域。

离散的傅里叶变换（DFT）

傅里叶变换：$F(u,v)={\sum }_{x=0}^{M-1}{\sum }_{y=0}^{N-1}f(x,y)e^{-j2\pi (ux/M+vy/N)}$
傅里叶反变换：$f(x,y)={\sum }_{u=0}^{M-1}{\sum }_{v=0}^{N-1}F(u,v)e^{j2\pi (ux/M+vy/N)}$
$x=0,1,2,...,M-1,y=0,1,2,...,N-1$

离散频率域中的相应间隔$\Delta u=\frac{1}{M\Delta T}和\Delta v=\frac{1}{N\Delta Z}$
频率域样本间的间隔与空间样本间的间距和样本数成反比。

（三）二维离散傅里叶变换的一些性质

• 平移性:傅里叶变换的平移性质表明了函数与一个指数相乘等于将变换后的频域中心或者空域中心移到新的位置。
  还可知，平移将不改变频谱的幅值。可以通过平移，将图像的原点移到屏幕中央。
• 周期性和共轭对称性：只需要一个周期里的变换就可以将F(u,v)在频率域完全确定。
  频谱是关于原点对称的
• 可分离性：可以通过先沿输入图像的每一列计算一维变换，然后沿中间结果的每一行再计算一维变换的方法来求二维变换。亦可先沿行再沿列进行。逆变换也可分离。
• 旋转性：$f(r,\theta +{\theta }_0)⟺F(w,\varphi +{\theta }_0)$
• 线性性：$F(a{f}_1(x,y)+b{f}_2(x,y))=aF(f_1(x,y))+bF(f_2(x,y))$
• 相似性（尺度变换）：$af(x,y)⟺aF(u,v)$
  $f(ax,by)⟺\frac{1}{|ab|}F(\frac{u}{a},\frac{v}{b})$

通常，相位角图像的视觉分析很少能得到直觉的信息，DFT的谱的分量决定正弦波的幅度，在任何给定的频率处，较大的幅度意味着图像中该频率的正弦波比较突出，反之，较小的幅度意味着图像中出现的正弦波较少。
$仅用相位角重建的图像$灰度信息会丢失，但图像中关键的形状特性是不会错的；
$仅用谱重建的图像$，仅包含灰度信息，直流项占支配地位，但图像中没有形状信息。

二、频率域滤波器

频率域中的滤波技术是以如下处理为基础的：修改傅里叶变换以达到特殊目的，然后计算IDFT返回到图像域。
基本的滤波公式为：$$g(x,y)={\zeta }^{-1}[H(u,v)F(u,v)]$$
其中${\zeta }^{-1}$是IDFT，$F(u,v)$是输入图像$f(x,y)$的DFT，$H(u,v)$是滤波函数，$g(x,y)$是滤波后的（输出）图像。

我们将使用低通滤波器模糊一幅图像，使用高通滤波器增强尖锐的细节。

（一）低通滤波器

$D_0$是截止频率
理想低通滤波器（ILPF）
$$H(u,v)=\{\begin{array}{ll}1,& \text{D(u,v)≤D0)}\\ 0,& \text{D(u,v)>D0)}\end{array}$$
巴特沃斯低通滤波器(BLPF)
$$H(u,v)=\frac{1}{1+[D(u,v)/D_0{]}^{2n}}$$
高斯低通滤波器（GLPF)
$$H(u,v)=e^{-D^2(u,v)/2{\sigma }^2},\sigma =D_0$$

（二）高通滤波器

理想高通滤波器（IHPF）
$$H(u,v)=\{\begin{array}{ll}0,& \text{D(u,v)≤D0)}\\ 1,& \text{D(u,v)>D0)}\end{array}$$
巴特沃斯高通滤波器(BHPF)
$$H(u,v)=\frac{1}{1+[D_0/D(u,v){]}^{2n}}$$
高斯高通滤波器（GHPF)
$$H(u,v)=1-e^{-D^2(u,v)/2{\sigma }^2},\sigma =D_0$$

（三）同态滤波

同态滤波器实现过程：

同态系统：通过非线性变换将非线性组合信号变换为线性组合，便于进行线性 处理的一类系统的总称。
同态滤波器是同态系统的一个特例。

（四）选择性滤波

常用的图像恢复方法有带阻滤波器、带通滤波器 、陷波滤波器等。
带阻滤波器：带阻滤波器消除或衰减了傅里叶变换原点附近的频段。
带通滤波器：带通滤波器执行与带阻滤波器相反的操作，可用全通滤波器减去带阻滤波器来实现带通滤波器。
陷波滤波器：陷波滤波器拒绝（或通过）事先定义的关于频率矩形中心的一个邻域的频率，可以选择性的修改DFT的局部区域

实践部分

原图：

截止频率为15

低通：

高通：

截止频率为30

低通：

高通：

带通滤波

总结(附代码)

1. 随着$D_0$的增大，低通滤波器的模糊水平降低，高通滤波器的锐化水平降低。
2. ILPF会产生振铃现象；与ILPF相比, BLPF变换传递函数比较平滑；GLPF图像边缘的模糊程度较BLPF产生的要大一些
3. IHPF通常具有振铃性质，微小物体和线条被加强，呈白色；BGPF比IHPF要平滑，边缘失真小；GHPF效果比BHPF稍弱，振铃现象不明显
4. 一阶BLPF没有振铃，二阶稍有但不明显
5. 高通滤波器让高频成分通过，使低频成分削弱，再经傅立叶逆变换得到边缘锐化的图像。$H_{hp}(u,v)=1-H_{lp}(u,v)$
6. 带通滤波使得图像的部分细节被突出

import cv2                #高通、低通滤波
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread("C:/Users/13121/Desktop/1.png", 0)
img=cv2.resize(img,(224,224),interpolation=cv2.INTER_CUBIC)#将原图缩小，方便观察
rows, cols = img.shape[:2]
nrows = cv2.getOptimalDFTSize(rows)# 计算最优尺寸
ncols = cv2.getOptimalDFTSize(cols)
nimg = np.zeros((nrows, ncols))# 根据新尺寸，建立新变换图像
nimg[:rows,:cols] = img
fft_mat = cv2.dft(np.float32(nimg), flags = cv2.DFT_COMPLEX_OUTPUT)#实部和虚部
fft_mat = np.fft.fftshift(fft_mat)# 反换位，低频部分移到中间，高频部分移到四周
d0 = 30

def main():
    imshow()
def fft_distances(m, n):# 计算D(u,v)
    u = np.array([i - m/2 for i in range(m)], dtype=np.float32)
    v = np.array([i - n/2 for i in range(n)], dtype=np.float32)
    ret = np.ones((m, n))
    for i in range(m):
        for j in range(n):
            ret[i][j] = np.sqrt(u[i]*u[i] + v[j]*v[j])
    u = np.array([i if i<=m/2 else m-i for i in range(m)], dtype=np.float32)
    v = np.array([i if i<=m/2 else m-i for i in range(m)], dtype=np.float32)
    return ret

def change_filter(flag):# 选择低通滤波器
    if flag == 1:    # 理想低通滤波器
        filter_mat = np.zeros((nrows, ncols ,2), np.float32)
        cv2.circle(filter_mat, (np.int(ncols/2), np.int(nrows/2)) , d0, (1,1,1), -1)        # 可以以(ncols/2, nrows/2)为圆心、d0为半径的圆内的值设置为1
    elif flag == 2:    # 巴特沃斯低通滤波
        n = 2 # 2阶
        filter_mat = None
        duv = fft_distances(*fft_mat.shape[:2])
        filter_mat = 1 / (1+ np.power(duv/d0, 2*n))
        filter_mat = cv2.merge((filter_mat, filter_mat))
    elif flag == 3:    # 高斯低通滤波（σ=d0）
        filter_mat = None
        duv = fft_distances(*fft_mat.shape[:2])
        filter_mat = np.exp(-(duv*duv) / (2*d0*d0))
        filter_mat = cv2.merge((filter_mat, filter_mat))
    return filter_mat
    
def ifft(fft_mat):# 对图像进行傅里叶反变换，返回反变换图像
    f_ishift_mat = np.fft.ifftshift(fft_mat)
    img_back = cv2.idft(f_ishift_mat)    # 傅里叶反变换
    img_back = cv2.magnitude(*cv2.split(img_back))    # 将复数转换为幅度, sqrt(re^2 + im^2)
    cv2.normalize(img_back, img_back, 0, 255, cv2.NORM_MINMAX)
    return np.uint8(np.around(img_back))[:rows,:cols]

def imshow():
    img1 = ifft(change_filter(1) * fft_mat) # 理想低通
    img2 = ifft(change_filter(2) * fft_mat) # 巴特沃斯低通
    img3 = ifft(change_filter(3) * fft_mat) # 高斯低通
    img4 = ifft((1 - change_filter(1)) * fft_mat) # 理想高通
    img5 = ifft((1 - change_filter(2)) * fft_mat) # 巴特沃斯高通
    img6 = ifft((1 - change_filter(3)) * fft_mat) # 高斯高通

    cv2.imshow('img.jpg',img)
    cv2.imshow('ILPF.jpg', img1)
    cv2.imshow('BLPF.jpg', img2)
    cv2.imshow('GLPF.jpg', img3)
    cv2.imshow('IHPF.jpg', img4)
    cv2.imshow('BHPF.jpg', img5)
    cv2.imshow('GHPF.jpg', img6)
    cv2.waitKey()
    cv2.destroyAllWindows()
if __name__ == "__main__":
     main()

import cv2		#带通滤波
import numpy as np
import math
img = cv2.imread("C:/Users/13121/Desktop/c.jpg", 0)
img = cv2.resize(img, (224, 224), interpolation=cv2.INTER_CUBIC)  # 将原图缩小，方便观察
fft_mat = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)  # 实部和虚部
fft_mat = np.fft.fftshift(fft_mat)  # 反换位，低频部分移到中间，高频部分移到四周
rows, cols = img.shape
nrows, ncols = int(rows / 2), int(cols / 2)
w = 20
d0 = 10
mask = np.ones((rows, cols, 2), np.uint8)
for i in range(0, rows):
    for j in range(0, cols):
        d = math.sqrt(pow(i - nrows, 2) + pow(j - ncols, 2))
        if d0 - w / 2 < d < d0 + w / 2:
            mask[i, j, 0] = mask[i, j, 1] = 0
        else:
            mask[i, j, 0] = mask[i, j, 1] = 1
f = fft_mat * mask# 掩膜图像和频谱图像乘积
fft_mat = np.fft.ifftshift(f)
iimg = cv2.idft(fft_mat)# 傅里叶逆变换
res = cv2.magnitude(iimg[:, :, 0], iimg[:, :, 1])

cv2.imshow('original', img)
cv2.imshow('lvbo', res)
cv2.waitKey()
cv2.destroyAllWindows()

参考代码1
参考代码2