数字图像学笔记 —— 17. 图像退化与复原（自适应滤波之「最小二乘方滤波」）

维纳滤波的缺点

维纳滤波（Wiener Filter），虽然是一种非常强大的退化图像还原算法，但是从实验过程我们也发现它存在着致命的缺陷，那就是要求输入退化系统的 $F(u,v)$ 是已知的。分析维纳滤波的理论，我们发现它是通过给未知的退化系统 $H(u,v)$ 输入原始数据 $F(u,v)$ 后得到退化后的数据 $G(u,v)$，然后通过梯度下降算法与约束条件 $MSE=(\hat{F}-F{)}^2$ 找出最接近 $H(u,v)$ 的近似解 $\hat{H}(u,v)$。

所以，在求解 $\hat{H}(u,v)$ 的过程中，它要求我们必须要预先知道输入退化系统 $H(u,v)$ 前的图像 $F(u,v)$ 的信息。但是我们在实际生活中，不一定都能拥有这样的条件，可以让我们顺利的找出 $\hat{H}(u,v)$，比方说对老照片的恢复。

拍摄时间未详的居里夫人照，图片来自网络。

对于上述图片来说，存在很明显的退化，而且我们也没有办法回到过去拍一张高清的照片，所以要想让照片得到修复，获得较好的视觉感受，维纳滤波在这里就无法使用了。

所以后续的研究者提出了一种仅依赖均值与方差便可以估算出最好复原效果的 「约束最小二乘方滤波（Constrained Least Squares Filtering）」。接下来我们来看看它是如何展现「化腐朽为神奇」的力量。

约束最小二乘方滤波

先从退化模型的一般形式（频率空间）出发，我们有如下公式：

$$G(u,v)=H(u,v)F(u,v)+N(u,v)$$

这是因为噪声函数 $D(u,v)$ 通常被认为与输入图像 $F(u,v)$ 相互独立，即 $N(u,v)$ 是与 $F(u,v)$ 不相关的噪声函数。因此，这个模型可以被看作是信号 $F(u,v)$ 通过退化系统 $H(u,v)$ 得到观测结果 $G(u,v)$ 后再加上噪声 $N(u,v)$。

接下来，我们可以使用最小二乘法来估计未知的退化函数 $H(u,v)$。最小二乘法是一种优化方法，可以通过最小化误差平方和来得到最优解。具体来说，对于给定的观测结果 $G(u,v)$ 和输入图像 $F(u,v)$，我们可以计算出它们之间的误差平方和：

$$E=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|G(u,v{)}^2-H(u,v)F(u,v){|}^2$$

我们的目标是找到一个 $H(u,v)$，使得误差平方和 $E$ 最小化。因此，我们可以通过求解下面的优化问题来得到最优的 $H(u,v)$

$$\underset{H(u,v)}{\min }E=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|G(u,v{)}^2-H(u,v)F(u,v){|}^2$$

然而，我们需要加上一些约束条件，以防止退化函数 $H(u,v)$ 变得过于复杂，从而导致过度拟合（overfitting）。一个常见的约束条件是 $H(u,v)$ 的均值为 $1$，方差为 ${\sigma }_H^2$，即：

$$\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}H(u,v)=1$$

$$\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|H(u,v)-\frac{1}{MN}\sum _{i=0}^{M-1}\sum _{j=0}^{N-1}H(i,j){|}^2={\sigma }_H^2$$

这个约束条件的意义是，我们期望退化函数 $H(u,v)$ 的均值为 $1$，这是因为退化函数的总能量应该保持不变。同时，方差 ${\sigma }_H^2$ 用于限制 $H(u,v)$ 的复杂度，从而避免过度拟合。

根据上述约束条件，我们可以得到一个带有约束的最小二乘优化问题：

$$\underset{H(u,v)}{\min }E=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|G(u,v)-H(u,v)F(u,v){|}^2$$

$$s.t.\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}H(u,v)=1$$

$$\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|H(u,v)-\frac{1}{MN}\sum _{i=0}^{M-1}\sum _{j=0}^{N-1}H(i,j){|}^2={\sigma }_H^2$$

上面这个优化问题可以通过拉格朗日乘子法来求解。我们首先构造一个拉格朗日函数：

$$L(H,{\lambda }_1,{\lambda }_2)=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|G(u,v)-H(u,v)F(u,v){|}^2+{\lambda }_1(\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}H(u,v)-1)+{\lambda }_2(\frac{1}{MN}\sum _{u_0}^{M-1}\sum _{v=0}^{N-1}|H(u,v)-\frac{1}{MN}\sum _{i=0}^{M-1}\sum _{j=0}^{N-1}H(i,j){|}^2-{\sigma }_H^2)$$

其中 ${\lambda }_1$ 和 ${\lambda }_2$ 是拉格朗日乘子。接下来，我们对 $L(H,{\lambda }_1,{\lambda }_2)$ 分别对 $H(u,v)$、${\lambda }_1$ 和 ${\lambda }_2$ 求偏导数，并令它们等于 $0$，可以得到下面的一组方程：

$$\frac{\partial L}{\partial {\lambda }_1}=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}H(u,v)-1=0$$

$$\frac{\partial L}{\partial H(u,v)}=-2F(u,v)(G(u,v)-H(u,v)F(u,v))+{\lambda }_1+2{\lambda }_2(H(u,v)-\frac{1}{MN}\sum _{i=0}^{M-1}\sum _{j=0}^{N-1}H(i,j))=0$$

$$\frac{\partial L}{\partial {\lambda }_2}=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|H(u,v)-\frac{1}{MN}\sum _{i=0}^{M-1}\sum _{j=0}^{N-1}H(i,j){|}^2-{\sigma }_H^2=0$$

解这个方程组可以得到最优的 $H(u,v)$。具体地，我们可以先把 ${\lambda }_1$ 和 ${\lambda }_2$ 消去，然后将 $H(u,v)$ 写成矩阵形式 $H=[h_1,h_2,...,h_L]$，其中 $L=MN$。接着，将偏导数为 $0$ 的方程写成矩阵形式 $A\text{h}=\text{b}$，其中 $A$ 和 $\text{b}$ 可以根据上述方程计算得到，$\text{h}$ 是一个列向量，它包含了所有的 $h_i$。最后，我们可以通过求解下面的矩阵方程来得到最优的 $\text{h}$：

$$(A^{T}A+\lambda I)\mathbf{h}=A^T\mathbf{g}$$

其中 $\text{g}$ 是 $G(u,v)$ 的列向量形式，$I$ 是单位矩阵，$\lambda$ 是一个正则化参数，用于平衡拟合误差和模型复杂度。

最终，我们可以得到最优的 $H(u,v)$：

$$H(u,v)=h_{uN+v}$$

其中 $h_i$ 是列向量 $\text{h}$ 中的第 $i$ 个元素。根据最优的 $H(u,v)$，我们可以得到修复后的图像 $F_{restored}(u,v)$：

Frestored(u,v)=F−1{G(u,v)/H(u,v)}F_{restored} (u, v) = \mathcal{F}^{-1} \{ G(u,v) / H(u,v) \} Frestored​(u,v)=F−1{G(u,v)/H(u,v)}

其中 ${\mathcal{F}}^{-1}$ 是傅里叶逆变换。由于在实际应用中，$H(u,v)$ 可能会变得过于复杂，从而导致过度拟合。为了避免这种情况，我们可以使用正则化技术，例如 Tikhonov 正则化，来控制模型的复杂度。具体来说，我们可以将优化问题改写为：

$$\underset{H(u,v)}{\min }E=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|G(u,v)-H(u,v)F(u,v){|}^2+\lambda \sum _{u=0}^{M-1}\sum _{v=0}^{N-1}|H(u,v){|}^2$$

其中 $\lambda$ 是一个正则化参数，用于控制 $H(u,v)$ 的平滑度。这样，我们可以得到下面的矩阵方程：

$$(A^{T}A+\lambda I)\text{h}=A^T\text{g}$$

解这个方程可以得到带有正则化项的最优 $H(u,v)$，从而得到修复后的图像 $F_{restored}(u,v)$。

总之，约束最小二乘方滤波是一种非常强大的图像修复技术，它可以通过仅依赖于均值和方差的约束条件来估计未知的退化函数，并从退化的图像中恢复出尽可能接近原始图像的图像。相比于维纳滤波，它更具有实用性和适用性。

给一个实际例子吧

上述内容实在太多，但是OpenCV已经帮我们减少了很多工作量，所以就可以简化为几个关键函数，所以我们就可以这样去调用它们

def ConstrainedLeastSquaresFiltering(image_path):# Load the image as grayscaleimg = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)# Add Gaussian noise to the imagenoisy_img = img + np.random.normal(0, 20, size=img.shape)# Define the degradation matrixdegradation_matrix = np.array([[0.9, 0.1, 0], [0.1, 0.8, 0.1], [0, 0.1, 0.9]])# Define the constraint matrixconstraint_matrix = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])# Define the regularization parameteralpha = 0.1# Compute the inverse of the degradation matrixdegradation_matrix_inv = np.linalg.inv(degradation_matrix)# Compute the constrained least squares estimatecls_estimate = cv2.filter2D(noisy_img, -1, degradation_matrix_inv)cls_estimate = cls_estimate - alpha * cv2.filter2D(cls_estimate, -1, constraint_matrix)# Convert the image to uint8 and clip it to [0, 255]cls_estimate = np.clip(cls_estimate, 0, 255).astype(np.uint8)# Show the original image, the noisy image, and the restored imagecv2.imshow("Original Image", img)cv2.imshow("Noisy Image", noisy_img)cv2.imshow("Restored Image", cls_estimate)cv2.waitKey(0)cv2.destroyAllWindows()

那么自然输出的效果就变成下面这个样子了

怎么样，是不是很厉害的技术呢？