复现Python代码请先安装PyTorch以及OpenCV

滤波器分类

常见的图像滤波器从是否为线性角度分为线性滤波器和非线性滤波器，非线性滤波器包括统计排序滤波器和形态学滤波器等。从空域还是频域角度又可以分为空域滤波和频域滤波。

本文仅涉及线性滤波器，并通过傅里叶变换和卷积定理说明空域滤波和频域滤波的等价性。

以下内容为脉络梳理以及代码实现，不包含证明等细节。

空域线性滤波

线性滤波器即对图像进行线性运算的滤波器。空域线性滤波器在图像$f(x,y)$上使用滤波器核(或窗口、模板等)$k$滑动执行加权求和操作，即卷积运算。其(离散)定义为：
$$(k\star f)(x,y)=\sum _{s=-a}^a\sum _{t=-b}^{b}k(s,t)f(x-s,y-t)$$

其直观过程如图(其他见conv_arithmetic)，下层为输入图像(扩充虚线部分使得输入输出尺寸相等)，移动的阴影为$3\times 3$卷积核，上层为输出图片：

例如，使用$3\times 3$的卷积核对图像进行均值滤波，卷积核为：
$$\frac{1}{9}\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right]$$

import torch
import torch.nn.functional as F
import cv2
import torchvision

x = cv2.cvtColor(cv2.imread("water_ouzel_n01601694.jpeg"), cv2.COLOR_BGR2GRAY)

x = (torch.from_numpy(x) / 255.0).unsqueeze_(0).unsqueeze_(0)

kernel = torch.tensor([[[[1, 1, 1], [1, 1, 1], [1, 1, 1]]]]) / 9.0

# pad:填充卷积核半径大小使得输出图像尺寸不变
filtered_x = F.conv2d(F.pad(x, (1, 1, 1, 1)), kernel)

torchvision.utils.save_image(x, "x.png")
torchvision.utils.save_image(filtered_x, "filtered_x.png")

频域线性滤波

根据傅里叶的思想，任何周期函数都可以表示为不同频率的三角函数的线性组合。

图像作为一种信号，也可以使用傅里叶变换从空域转换到频率域。空域和频域的表示是等价且可以相互转换的。
二维傅里叶变换及其逆变换为：
$$\begin{gathered} F(u,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(t,z)e^{-j2\pi (ut+vz)}dtdz \\ f(t,z)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }F(u,v)e^{j2\pi (ut+vz)}dudv \end{gathered}$$

二维离散傅里叶变换及其逆变换为：
$$\begin{gathered} 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)=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{j2\pi (ux/M+vy/N)} \end{gathered}$$

使用欧拉公式
$$e^{j\theta }=\cos \theta +j\sin \theta$$
可以得到傅里叶变换结果在复数域下的表示：
$$F(u,v)=R(u,v)+jI(u,v)=|F(u,v)|e^{j\varphi (u,v)}$$

其中，$|F(u,v)|=\sqrt{R^2(u,v)+I^2(u,v)}$称为傅里叶频谱，$\varphi (u,v)=\arctan \left[\frac{I(u,v)}{R(u,v)}\right]$为相角或相位谱。

对左图进行傅里叶变换得到频谱图，以频谱图中心为原点，距离中心点越远的位置代表频率越低，其中：$F(0,0)={\sum }_{x=0}^{M-1}{\sum }_{y=0}^{N-1}f(x,y)=MN\bar{f}$，表示$f(x,y)$图像的均值：

注意：

1. 频谱图和原图的对应位置的像素没有任何关系，频谱图的像素位置代表其频率，像素值代表频谱或相角。
2. 通常在可视化频谱图时会对其幅度取对数——log，这是因为低频分量和高频分量强度呈现幂级关系，相关搜索：image structure & power law。

对频谱图进行中心化，使其更容易观察，此时，低频信息集中在中心点：

进行高通滤波时， 只需要去掉中心的低频信息再变换回空间域，就实现了图像在频域的高通滤波。
去除中心点低频信息后的频谱图：

转换回空域，就得到了高通滤波的结果：

import torch
import cv2
import numpy as np

spaital_x = cv2.imread("water_ouzel_n01601694.jpeg")

# dft
fre_x = torch.fft.fft2(torch.from_numpy(spaital_x).permute(2, 0, 1))

# high-pass filter
shift_fre_x = torch.fft.fftshift(fre_x)
c, w, h = shift_fre_x.shape
shift_fre_x[:, w//2-32:w//2+32, w//2-32:w//2+32] = 0
fre_x = torch.fft.ifftshift(shift_fre_x)

# idft
ispaital_x = torch.fft.ifft2(fre_x)

cv2.imshow("original", spaital_x)
cv2.imshow("m", 20 * torch.log(torch.abs(fre_x)).permute(1, 2, 0).numpy().astype(np.uint8))
cv2.imshow("shifted m", 20 * torch.log(torch.abs(shift_fre_x)).permute(1, 2, 0).numpy().astype(np.uint8))
cv2.imshow("ioriginal", torch.abs(ispaital_x).permute(1, 2, 0).numpy().astype(np.uint8))

cv2.waitKey(0)
cv2.destroyAllWindows()

卷积定理

空域中两个函数的卷积的傅里叶变换等于频域中两个函数的傅里叶变换的乘积。即空间卷积操作可以用频域上的乘法实现，这可以对卷积操作进行加速：
$$(f\star h)(x,y)\iff (F\cdot H)(u,v)$$

由于离散傅里叶变换本身固有的周期性，为了避免由此产生的交叠错误，需要在两个输入A长度和B长度样本(比如输入图像和卷积核)中填充零，使得他们长度相同P，且：
$$P\ge A+B-1$$

使用DFT实现卷积操作：

import torch
import cv2
import torch.nn.functional as F
import torchvision

# 自定义单通道卷积运算
def conv2d_fft(x: torch.Tensor, kernel: torch.Tensor):
    kw, kh = kernel.shape
    
    print(kernel.shape)
    
    kernel = kernel.rot90().rot90()
    
    # 填充零，避免周期性产生的交叠错误
    x = F.pad(x, [0, kw - 1, 0, kh - 1])
    k = F.pad(kernel, [0, x.size(0) - kw, 0, x.size(1) - kh])
    
    # DFT
    fre_x = torch.fft.fft2(x)
    fre_k = torch.fft.fft2(k)
    
    # IDFT
    return torch.real(torch.fft.ifft2(fre_x * fre_k))[kw//2:-(kw//2), kh//2:-(kh//2)].contiguous()


spaital_x = cv2.imread("water_ouzel_n01601694.jpeg")
spaital_x = torch.from_numpy(cv2.cvtColor(spaital_x, cv2.COLOR_BGR2GRAY)) / 255.0

kernel = torch.tensor([
    [-1, 0, 1], 
    [-2, 0, 2], 
    [-1, 0, 1]
], dtype=torch.float32)

y1 = conv2d_fft(spaital_x, kernel)

pad = (kernel.size(-1) // 2, kernel.size(-1) // 2, kernel.size(-2) // 2, kernel.size(-2) // 2)
y2 = F.conv2d(F.pad(spaital_x.view(1, 1, *spaital_x.shape), pad), kernel.view(1, 1, *kernel.shape))

print(y1)
print(y2)

torchvision.utils.save_image(torch.abs(y1), 'y1.png', normalize=True)
torchvision.utils.save_image(torch.abs(y2), 'y2.png', normalize=True)

空域滤波和频域滤波的关系

由于空间域滤波和频率域滤波之间可以通过卷积定理等价，那么也就意味着，空间域中的卷积核$h(x,y)$通过傅里叶变换就得到了频率域中对应等价的滤波器传递函数$H(u,v)$，相应的，已知$F(u,v)$时，也可以通过傅里叶反变换得到空间域中等效卷积核$h(x,y)$。两个滤波器之间构成一个傅里叶变换对：
$$h(x,y)\iff H(u,v)$$
由于在频率域中可以直观的处理频率，因此，在频率域设计合适的滤波器传递函数后直接滤波，或通过傅里叶反变换得到空间域等效核。

此外，快速傅里叶变换是实现卷积运算的一种非常高效方法

空域和频域的高斯滤波

频域中，高斯低通滤波器的传递函数为：
$$H(u,v)=e^{-\frac{D^2(u,v)}{2{\sigma }^2}}$$
其分辨率和输入图片一致，$\sigma =32$时如下图所示

与输入图片的乘积即为频率域的高斯低通滤波，保留了频域中更多的低频成分。将该传递函数使用傅里叶变换转换到空间域后，得到的卷积核如下：

频率域高斯函数的傅里叶反变换也是高斯函数且频率域中窄高斯传递函数意味着空间域中较宽的核函数。

空域和频域的边缘检测

空间域的边缘检测常用线性算子：一阶差分算子、Sobel算子、Prewitt算子、Laplace算子等：
$$\left[\begin{array}{ccc}0& 0& 0\\ -1& 0& 1\\ 0& 0& 0\end{array}\right]\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]\left[\begin{array}{ccc}-1& -1& -1\\ -1& 8& -1\\ -1& -1& -1\end{array}\right]$$
其转换到频率域后的频谱图分别为(分辨率:303 $\times$ 303):


均对低频进行了抑制，通过了高频频段。Sobel算子 效果好于一阶差分算子和Prewitt算子则是因为Sobel算子对图像各种方向的超高频信息进行了抑制，这对降低噪声的影响很有帮助。

此外，差分算子的间隔对带通频率也有影响：
$$\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& -1& 1& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& -1& 0& 1& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ -1& 0& 0& 0& 1\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ -2& -1& 0& 1& 2\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ -1& 0& 0& 0& 1\\ -2& -1& 0& 1& 2\\ -1& 0& 0& 0& 1\\ 0& 0& 0& 0& 0\end{array}\right]$$

为了抑制高频噪声，先对图像进行高斯滤波，然后再进行一阶导数，这等价于先对高斯求导再与图像进行卷积运算，二维高斯函数及对$x$方向的偏导数的如下：
$$\begin{array}{rl}G(x,y)& =\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{x^2+y^2}{2{\sigma }^2}}\\ \frac{\partial G}{\partial x}& =(-\frac{1}{2\pi {\sigma }^4})x{e}^{-\frac{x^2+y^2}{2{\sigma }^2}}\end{array}$$
方差为$1.0$的离散卷积核：
$$\left[\begin{array}{ccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -3& -8& 0& 8& 3& 0& 0\\ 0& -2& -16& -36& 0& 36& 16& 2& 0\\ 0& -3& -27& -60& 0& 60& 27& 3& 0\\ 0& -2& -16& -36& 0& 36& 16& 2& 0\\ 0& 0& -3& -8& 0& 8& 3& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$$
其对应频谱图为，明显的抑制了高频信号：